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SUMMARY 


The  need  that  led  to  this  research  project  was  to  obtain  mojre 
information  about  enqines  from  existing  data  without  inordinate  increase  in 
computation.  The  primary  focus  was  on  the  two  fundamental  parts  of  engine 
management,  the  estimation  of  engine  removal  times  and  the  use  of  these 
estimates  to  predict  replacement  requirements.  The  secondary  focus  was  on 
procedures  for  estimating  spare  engine  requirements,  improving  simulation, 
and  updating  estimators. 

The  actuarial  method  has  served  long  and  well  to  provide  fundamental 
information  for  engine  management  decisions,  and  technological  development 
has  progressed  to  the  point  where  improvement  in  the  accuracy  of  actuarial 
estimates  can  be  achieved  without  a major  increase  in  computation.  The 
statisticians  dream  of  more  information  from  the  same  data  can  be  achieved 
because  all  the  information  now  available  is  not  used.  For  the  sake  of 
computational  convenience,  information  collected  on  engine  removal  times  is 
summarized  in  counts  of  the  numbers  of  engine  exposures  and  removals  in 
actuarial  age  intervals.  Theoretical  developments  now  allow  the  engine 
removal  times  themselves  to  be  used  for  estimation  of  removal  times  and  for 
prediction  of  replacement  requirements  without  inordinate  computation. 

Many  recommendations  for  improvement  of  the  actuarial  method  and  for 
improvement  of  engine  management  follow  from  the  developments  that  will  be 
described  in  this  report.  Some  are  recommendations  for  development  of 
promising  improvements  and  some  are  changes  that  can  be  made  now  permanently 
or  perhaps  temporarily  until  further  improvement  is  possible.  The  nature 
and  immediacy  of  each  recommendation  will  be  indicated. 
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SECTION  I 

REVALUATION  OF  THE  ACTUARIAL  METHOD 
1.  Summary  of  Revaluation 

1 

t There  are  several  potential  improvements  that  could  be  made  in  estima- 

tion of  engine  removal  times  and  failure  rates.  Some  can  be  done  now, 
some  require  further  study,  testing  and  programming,  and  some  involve 
eventual  replacement  of  the  actuarial  failure  rate  estimator  with  an 
alternative  estimator.  1 

The  first  recommendation  is  not  to  smooth  the  crude  failure  rates  and 
not  to  throw  out  data  where  data  is  sparse.  The  latter  suggestion  raises 
the  question  of  scarce  data  and  the  confidence  in  very  low  failure  rate 
estimates.  This  question  can  be  resolved  by  recommendation  two,  an  actuarial 
type  estimator  that  uses  "flexible"  age  intervals  that  contain  either  equal 
numbers  of  exposures  or  removals  in  each  age  interval.  The  actuarial 
estimator  with  flexible  intervals  was  not  used  for  comparison  with  the 
currently  used  crude  actuarial  estimators  because  still  better  estimators 
were  available. 

The  actuarial  estimator  of  failure  rates  summarizes  the  data  on  engine 
removal  times,  operating  times,  and  ages  of  operating  engines  into  counts 
of  the  numbers  of  exposures  and  removals  in  each  actuarial  age  interval. 

This  summary  of  data  reduces  the  amount  of  information  in  the  sample  and 
consequently  the  accuracy  of  any  estimator.  An  estimator  that  uses  the 

. 

actual  engine  removal  times  should  be  better  than  an  estimator  that  uses 
only  counts  of  the  numbers  of  removals  during  age  intervals.  An  estimator 
that  also  incorporates  the  ages  of  engines  still  operating  should  be 
better  yet.  And  an  estimator  that  incorporates  the  operating  time  during 

** 

a period  should  give  the  sharpest  picture  of  engine  operation  during 

I 
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that  period.  An  estimator  that  uses  all  three  types  of  information  has  not 
been  derived  yet,  but  that  is  the  next  step  for  research. 

The  "product  limit"  estimator  of  Kaplan  and  Meier  [ 4 ] uses  the 
removal  times  and  the  ages  of  engines  still  in  operation  to  give  a dis- 
tribution free  estimate  of  the  cumulative  distribution  function  of  engine 
removal  times.  It  can  be  adapted  to  also  give  a cumulative  failure  rate  function 
estimator  which  in  turn  can  be  used  to  obtain  interval  failure  rates 
comparable  to  actuarial  failure  rates  for  any  desired  interval  widths. 

For  one  additional  assumption  on  the  underlying  nature  of  engine 
removals,  a maximum  likelihood  estimator  of  the  cumulative  distribution 
function  is  derived  in  Appendix  A that  is  simpler  than  the  product  limit 
estimator  and  can  also  be  used  in  engine  diagnosis.  The  assumption  is  that 
inspection  removals  and  usage  removals  are  statistically  independent,  a 
testable  hypothesis. 

All  three  estimators,  here  called  the  actuarial  estimator,  the  product 
limit  estimator,  and  the  maximum  likelihood  estimator  will  be  compared  on 
the  basis  of  how  closely  they  estimate  the  theoretical  distribution  of  simulated 

engine  removal  time  data  and  on  the  basis  of  how  well  they  predict  replace- 
ment requirements.  On  theoretical  grounds,  the  actuarial  estimator  is  the 
least  accurate  because  all  three  estimators  are  in  fact  maximum  likelihood, 

Barlow  and  Proschan  [ 5 ],  but  for  different  sample  information  and  for 
different  engine  removal  time  models.  The  actuarial  estimator  uses  summary  counts 
of  exposures  and  removals.  The  actuarial  and  product  limit  estimators  make 
no  assumption  about  the  engine  removal  time  mechanism.  The  product  limit 
and  maximum  likelihood  estimators  use  the  removal  time  data  and  the  ages 
of  surviving  engines,  but  the  maximum  likelihood  estimator  is  derived  from 
a removal  time  model  that  assumes  statistical  independence  between  inspection  and 
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usage  removal  times.  Whenever  a random  process  is  known  to  have  more 
statistical  structure,  more  information  can  be  inferred  from  data  about  the 
process,  so  it  is  believed  that  the  maximum  likelihood  estimator  is  the  most 
efficient  estimator.  Analytical  proof  of  this  has  not  been  completed. 

It  is  only  fair  to  point  that  the  existing  actuarial  method  has 
provision  for  counting  only  flying  experience  during  a limited  period  of 
time.  The  product  limit  and  maximum  likelihood  estimators  use  removal  time 
data  and  ages  on  engines  which  may  have  begun  operation  at  any  time  in  the 
past.  Comparison  of  all  three  estimators  was  made  on  an  equal  basis,  all 
engines  new  at  the  start  of  the  data  gathering  process  although  all  may  not 
have  failed  by  the  end  of  the  data  collection  period.  In  this  comparison, 
the  product  limit  and  maximum  likelihood  estimators  were  noticeably  better 
even  for  a very  small  sample.  However,  the  actuarial  estimator  will  probably 
yield  a better  answer  to  questions  of  comparing  the  experience  of  one 
calendar  quarter  to  another.  Research  will  be  continued  to  obtain  Droduct 
limit  and  maximum  likelihood  estimators  freely  comparable  to  the  actuarial 
estimator. 

An  updating  program  has  been  written  that  compares  current  data  on 
engine  removal  times  and  ages  of  engines  that  continue  to  operate  with 
past  data  of  the  same  type.  It  provides  a measure  of  how  likely  it  is  that 
the  current  and  all  past  data  subsets  came  from  populations  of  engine 
removal  data  that  are  the  same.  These  measures  will  provide  the  Aerospace 
Engine  Life  Committee  with  guidance  in  selection  of  data  to  be  included 
in  the  Official  Failure  Rate.  The  program  uses  a Wilcoxon  type  rank  test 
statistic  first  proposed  by  Geuan  [ 6 ].  This  program  is  proposed  as  an 
alternative  to  the  currently  used  "statistical  test"  for  determining  whether 


the  actuarial  failure  rate  during  an  age  interval  has  changed.  The 
program  is  subject  to  the  same  qualification  mentioned  earlier,  that  the 
actuarial  test  compares  only  experience  obtained  during  the  current  calendar 
quarter  while  the  updating  program  compares  a removal  time  data  set  that  may 
include  engines  which  had  begun  operation  prior  to  the  current  quarter. 

It  remains  to  be  seen  whether  this  distinction  is  important.  Meanwhile 
research  has  been  started  on  a modified  Gehan  type  updating  program  that 
compares  current  quarter  experience  with  past  experience.  Even  so,  the 
Wilcoxon  and  Gehan  type  rank  test  is  a temporary  alternative  proposed  for 
use  only  until  an  updating  program  based  on  the  product  limit  or  maximum 
likelihood  estimators  can  be  developed  and  until  either  the  question  of 
use  of  current  quarter  experience  can  be  resolved  or  until  the  product  limit 
and  maximum  likelihood  estimators  can  be  modified  to  incorporate  engine 
operation  experience  from  a fixed  calendar  interval. 
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2.  Review  of  the  Actuarial  Estimator  of  Failure  Rates 


The  actuarial  estimator  of  the  failure  rate  during  an  engine  ooerating 
age  interval  is  basically  the  number  of  engines  with  operating  age  at 
removal  recorded  durinq  that  interval  divided  by  the  number  of  engines  that 
operated  until  the  beginning  of  that  age  interval.  Let  At  be  the  width  of 
an  age  interval,  d-  the  number  of  engines  removed  during  the  interval 
with  ages  between  (i-1)  At  and  iAt,  and  n.  the  number  of  engines  that 
operated  at  least  until  age  (i-l)At,  i=l,2...,k,  then  the  actuarial 
failure  rate  estimator  is 

Ri  = V"i  i = 1 ,2, . . . ,k.  (2.1) 

These  simple  estimators  must  be  modified  to  account  for  data 
on  engines  that  have  been  operated  but  not  yet  failed  at  the  time  data  was 
collected  otherwise  the  actuarial  failure  rate  estimators  will  be  biased 
high.  This  additional  information  referred  to  here  as  survivors'  ages  is  used 
to  modify  the  denominator  n^  Dy  counting  the  total  numbers  of  engines, 
survivors  and  removed,  entering  each  age  interval  and  counting  a fractional 
exposure  in  the  age  interval  if  a survivor's  age  is  part  of  the  way  through 
an  age  interval.  E.g.  if  an  engines'  age  at  time  of  data  collection  was  750 
hours  and  the  actuarial  interval  was  from  720  to  760,  three  quarters  of 
an  exposure  would  be  credited  to  the  age  interval  and  number  of  exposures  n^ 
in  that  interval. 

A further  modification  is  used  if  it  is  desired  to  estimate  failure 
rates  based  only  on  flying  experience  during  a particular  calendar  period. 
Fractional  exposures  are  counted  for  engines  whose  age  at  the  beginning  of 
the  period  was  part  way  through  an  age  interval,  and  a full  exposure 
is  counted  for  engines  removed  during  an  interval  even  though  it  may  have 
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begun  the  calendar  period  with  operating  age  in  the  same  actuarial  age 
interval  within  which  it  was  removed.  The  resulting  ratio  of  deaths  to 
exposures  is  the  "crude  actuarial  failure  rate"  estimator  of  T.O.  00-25-128 
[7  ] and  T.O.  00-25-217  [ 8 ]. 

The  crude  actuarial  failure  rates  are  then  smoothed  with  a weighted 
moving  average  up  to  age  intervals  where  there  are  insufficient  exposures 
and/or  removals.  Thereafter,  actuarial  failure  rates  are  extrapolated  by 
linear  regression  to  the  maximum  operating  time. 

The  crude  failure  rates  (2.1)  have  several  desirable  large  sample 
statistical  properties  - 

1.  They  are  maximum  likelihood  estimators  for  the  sample 
information  { ( d^ , n^),  i = l,2,...,k}. 

2.  They  are  asymptotically  jointly  normally  distributed  with 
mean  zero  and,  surprisingly,  with  zero  covariance  as  the 
sample  size  grows  under  the  condition  that  the  probability 
of  an  engine  being  removed  at  age  zero  is  zero. 

The  first  assertion  follows  from  a related  assertion  in  Barlow  and 
Proschan  [ 5 ] which  has  been  known  for  many  years,  and  the  second  assertion 

follows  from  theorem  1.2  in  Barlow  and  Proschan  [ 5 ] which  is  proved  for 

an  estimator  that  is  R^/At. 

The  actuarial  estimator  was  developed  for  use  in  the  life  insurance 
industry  for  predicting  human  failure  rates.  A typical  human  failure  rate 
estimate  is  shown  in  figure  1. 


re  Rajtes  far  the  .Total  Population  ot  the  United  States  in  1P59- 19161 


It  is  inherently  a smooth  function  except  in  the  first  few  years  because 
there  is  nothing  in  the  mechanism  causing  deaths  that  is  age  specific  such 
as  components  that  wear  out  at  specific  ages.  Consequently,  actuaries 
dependent  on  sampling  to  construct  their  failure  rates  feel  justified 
in  smoothing  any  irregularities  that  appear  in  their  failure  rate  estimators. 
And  they  have  larger  samples  than  are  available  to  the  Air  Force  for  each 
engine  type.  The  weighted  moving  average  smoothing  technique  has  been 
borrowed  from  the  actuarial  industry  to  smooth  the  crude  engine  failure  rates 
for  the  following  explanation  (T.O.  00-25-123  p 4-1  [ 7 ]). 

"As  estimated,  these  rates  are  subject  to  statistical  error  which 
may  be  large  or  small,  depending  upon  chance  and  the  volume  of  the 
data  from  which  the  crude  rates  were  computed.  The  error  is  as 
likely  to  he  positive  as  negative , i.e.,  the  probability  of  the 
estimate  being  high  is  equal  to  the  probability  of  the  estimate  being 
low.  By  inter-relating  the  rates  for  a number  of  adjacent  age  intervals 
the  error  for  each  rate  is  redistributed  among  all  the  rates  in  the 
group.  When  this  is  done  the  positive  errors  tend  to  balance  with 
negative  errors,  and  all  unusually  high  errors  are  spread  over  all 
of  the  rates  so  that  the  crude  rates  are  brought  closer  to  their  true 
rates. " 

The  weights  in  the  moving  average  smoothing  formula  vary  depending  on  the 
"failure  density",  the  average  number  of  removals  per  interval  in  the  range 
where  smoothing  is  to  be  applied  (T.O.  00-25-217  p.  7-1 1 [ 8 ]).  A 7- 
point  smoothing  formula  is 
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R.  = (-.0587)  (R._3+  Ri+3)  + 
+ .0587  (R. _2+  Ri+2)  + 
+ .2937  (Ri _-|  + R.+1)  + 


+ .4126  R. 


(2.2) 


It  is  recommended  that  smoothing  be  stopped  immediately  for  the 
following  reasons: 

1.  The  underlying  mechanism  of  engine  removals  contains  factors 
that  create  major  changes  in  failure  rates  from  one  interval 

to  the  next,  factors  such  as  the  high  probability  of  removal  at 
inspection  and  the  possibility  that  certain  engine  components  may 
fail  within  an  age  interval.  Smoothing  suppresses  the  infor- 
mation contained  in  inspection  removal  times  and  useful 
diagnostic  information  of  high  usage  failure  rates  possibly 
attributable  to  particular  engine  components. 

2.  Smoothing  is  appropriate  primarily  to  time  series  such  as  stock 
market  prices  or  human  ages  which  advance  on  or  close  to  a one 
for  one  relation  with  calendar  time.  Engine  ages  do  not 
advance  in  direct  proportion  to  calendar  time  and  so  smoothing  an 
interval  failure  rate  over  neighboring  intervals  is  not  appropriate. 

3.  The  explanation  (italics)  that  errors  are  as  likely  to  be  high  as 
low  is  not  known  by  the  author  to  be  true  except  asymptotically. 

4.  The  explanation  (italics)  that  positive  errors  tend  to  cancel 
negative  errors  is  fallacious.  Due  to  reason  1,  it  is  not 
known  whether  a high  crude  failure  rate  is  in  error  to  the  high 
side.  The  summarization  of  removal  time.,  into  counts  of  removals 
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per  age  interval  may  mean  that  a high  failure  rate  is  in 
error  on  the  low  side  of  the  true  value! 

5.  The  negative  weighting  factors  in  the  smoothing  formula 
(2.2)  are  appropriate  only  when  there  is  a detected  cyclic 
error  in  the  failure  rates,  i.e.,  it  is  known  that  the 
crude  failure  rates  3 age  intervals  ago  and  3 intervals  ahead 
(for  the  7-point  smoothing  formula)  are  in  error  in  the 
opposite  direction  to  the  central  age  interval's  failure  rate. 

This  is  almost  assuredly  not  true  because  the  asymptotic  covariance 
of  failure  rates  is  zero  (Barlow  and  Proschan  [ 5 ])  and  because 
there  is  no  cyclic  mechanism  in  the  engine  failure  process  with 
such  periodic  behavior.  Furthermore,  rather  sophisticated  means 
are  required  to  detect  cyclic  behavior  and  do  not  call  for  a 
changing  period  for  the  negative  weight  as  occurs  when  changing 
from  the  7-point  to  another  smoothing  formula. 

Despite  this  condemnation  of  the  smoothing  procedure,  the  crude  failure 
rates  are  reasonable  and  convenient  estimators,  and  are  easier  to  compute 
than  smooth  failure  rates.  Objection  number  1 regarding  the  underlying 
mechanism  of  engine  failure  has  led  to  the  proposal  of  an  engine  removal 
time  model,  section  4,  for  which  the  maximum  likelihood  estimator  of  the 
usage  removal  time  cumulative  distribution  function  and  inspection  removal 
probabilities  have  been  derived,  Appendix  A.  That  estimator  is  based  on 
the  actual  removal  times  and  survivor's  ages.  Its  performance  is  compared 
to  the  crude  actuarial  failure  rate  estimator  in  section  5,  and  it  is  found 
to  be  superior.  Nevertheless,  the  crude  failure  rate  estimator  may  be 
improved  by  changing  to  flexible  age  intervals  as  described  in  the  next  section. 
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3.  Variable  Actuarial  Age  Intervals 


Examination  of  the  variance  of  the  percent  survivor  function  or  cumula- 
tive distribution  function,  resulted  in  the  conclusion  that  its  variance 
could  be  minimized  by  having  the  same  number  of  exposures  in  each  age  inter- 
val! The  variance  of  any  estimate  is  a measure  of  the  precision  of  that  estimate, 
and,  if  the  estimate  is  unbiased,  then  the  estimate  with  lowest  variance  is 
also  the  most  accurate. 

The  conclusion  about  minimizing  variance  as  a function  of  the  number 
of  exposures  was  arrived  at  indirectly  by  examining  the  variance  of  an 
actuarial  type  estimate  of  the  percent  survivor  function  or  tail  cumulative 
distribution  function.  First,  some  notation  must  be  defined,  and  the  condi- 
tions for  the  derivation  must  be  outlined.  The  definitions  and  formulas 
are  taken  from  Gross  and  Clark  [3  ] Chapter  2. 

{t.j}  is  the  sequence  of  age  interval  boundaries  i =0,1 k+l  where 

tk+1  may  be  », 

th 

{dj}  is  the  sequence  of  numbers  of  removals  in  the  iLn  interval  ( t^ _-j , 
t,), 

{n-j}  is  the  sequence  of  numbers  of  engine  exposures  during  the  i^ 
interval,  and 

{ P ( t , ) > is  the  sequence  of  estimated  percent  surviving  to  time  t^  esti- 
mated as  the  cumulative  proportion  in  the  sample  data 

P(t0)  = 1.0,  P(t| ) = etC-  P(ti}  = di-l/ni-T^ti-l) 

(The  true  unknown  value  of  percent  surviving  is  denoted  P(ti)). 

It  should  be  noted  that  the  authors  Gross  and  Clark  assume  that  partial 
engine  exposures  are  uniformly  distributed  over  each  interval  so  that  the 
number  of  exposures  is  (approximately)  equal  to  the  number  of  engines  that 
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enter  an  age  interval  minus  half  the  number  of  engines  that  have  partial 
exposures  during  the  interval  (other  than  engines  removed  during  the  inter- 
val that  are  counted  as  a full  exposures).  Consequently,  conclusions  re- 
garding the  variance  of  estimators  as  a function  of  the  number  of  exposures 
per  interval  are  indicative  rather  than  absolutely  accurate  conclusions. 

In  fact,  the  variance  formula  to  be  used  in  obtaining  the  conclusion  is 
itself  an  approximation. 

A 

The  variance  of  P( t ^ ) , the  estimate  of  the  percent  surviving  to  t^, 
is  approximately  (Gross  and  Clark  p.  41  [3]) 


i-1 


Var  P(ti)  s [P(ti )]  ^ dj/(nj(nj-dj)  ) 


(1) 


as  a function  of  both  the  sequences  {n j } and  {d j } . The  numbers  of  removals 
d,  in  each  interval  will  be  assumed  equal  to  1.0  for  convenience  in  showing 
that  the  variance  can  be  minimized  by  choosing  intervals  to  contain  equal 
numbers  of  exposures  n^.  (It  is  recognized  that  these  are  possibly  incom- 

J 

patible  events  unlikely  to  occur  except  by  accident,  but  the  conclusion  is 
used  only  as  an  indication).  Then  the  problem  is  to  minimize  the  variance 

A 

of  P(t-j)  by  choice  of  the  numbers  of  exposures  to  be  included  in  each  age 
interval  (by  adjusting  the  age  interval  widths). 


min  = Var[P(ti)]  = min  [P(ti)]2  l (n jKJ)"1 

i^l  > n^»  . • . j n i_i  j ~ 1 


where  the  {n.}  are  constrained  by  the  requirement  that 
J 

i-1 

I n,  = N(t-) 
j=l  J 

where  N(t^)  is  the  number  of  exposures  up  to  time  t^. 
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Constrained  minimization  can  be  done  by  the  LaGrange  multiplier  method 
by  replacing  the  variance  minimization  by  the  LaGrangian.  ( [P ( t ^ ) ] has 
been  omitted  from  the  minimization  because  it  is  constant). 

i-1  , i-1 

min  l (n,(n  • + X ( J n.-N(ti ) ) = min  L 

rip  n2  ....  ni  p A j=l  J J j=l 


where  A is  the  "LaGrange  multiplier".  Calculus  is  used  to  find  the  minimum, 
|j-  = -(taj-D/Injlnj.,)  )2  * X - 0,  j-1,  2 i-1 

I r-j’  -J-NCt,)  - 0 

Since  all  the  derivatives  6L/6n j yield  the  same  equation  subject  to  the 
constraint 

i-1 

l n-  = N(tJ, 
j=l  J 

the  optimal  allocation  of  exposures  to  intervals  is  to  choose  equal  numbers 
n^  in  each  interval  (as  closely  as  possible  because  N(t-j)  may  not  be  evenly 
divisible  by  i-1 ) . 

The  conclusion  is  that  to  reduce  the  variance  of  the  estimate  of  a 
particular  percent  surviving,  choose  age  intervals  to  contain  the  same  number 
of  exposures.  This  indicates  that  the  age  intervals  should  be  short  where 
there  are  many  exposures  and  long  where  there  are  few.  Since  many  exposures 
usually  result  in  more  removals,  the  conclusion  above  also  implies  that  age 
intervals  should  be  short  where  there  are  many  removals  and  long  where  there 
are  few.  In  fact,  if  the  installation  and  operation  of  engines  makes  the 
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number  of  exposures  sequence  {n.}  mathematically  independent  of  the  sequence 

J 

A 

of  removals  {d  . } then  the  minimization  of  Var  P(t.j)  with  respect  to  the 
sequence  (d j } yields  the  same  conclusion,  that  the  age  intervals  should  be 
chosen  to  have  an  equal  number  of  removals  in  each  interval  to  minimize 
Var  Pftj). 

The  "product  limit"  estimator  of  Kaplan  and  Meier  [4]  carries  the 
concept  of  equal  numbers  of  removals  per  interval  to  its  ultimate  extension, 
one  removal  per  interval.  Further  investigation  of  the  optimal  interval 
widths  to  minimize  variance  was  not  done  because  this  ultimate  extension  was 
available  and  because  the  data  on  engine  removal  times  and  ages  was  already 
at  hand.  However,  the  calculus  method  used  in  optimization  above  may  be 
used  to  choose  interval  widths  for  actuarial  estimators  if  it  is  desired  to 
continue  to  use  actuarial  estimators.  Further  work  is  required  however 
since  a computer  program  has  not  been  written  for  actuarial  estimators  with 
equal  numbers  of  removals  or  exposures  per  age  interval  and  since  the 
simplifying  assumptions  made  in  the  derivation  may  not  be  true  in  practice. 
Additional  calculation  is  needed  to  obtain  optimal  age  intervals  when  removals 
and  exposures  are  mathematically  dependent. 
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4.  Alternative  Estimators 


Examination  of  crude  failure  rates  for  many  types  of  engines  shows  high 
failure  rates  in  the  age  intervals  surrounding  inspection  times.  This  is 

because  the  elaborate  testing  done  at  inspection  is  much  more  likely  to 
disclose  an  impending  cause  for  removal  than  during  usage.  On  the  other  hand, 
there  is  reason  to  believe  that  normal  maintenance  carried  out  at  inspections 
serves  only  to  permit  an  engine  to  achieve  its  normal  operating  age  at 
eventual  removal  for  repair  or  rebuilding.  Inspection  maintenance  seems  no 
more  likely  to  rejuvenate  an  engine  than  to  damage  it  further.  These  obser- 
vations led  to  proposal  of  a statistical  model  of  engine  removal  times  that 
assumes  independence  of  inspection  and  usage  removal  times  (George  [9]). 

This  model  assumes  a cumulative  distribution  function  F^(t)  for  usage  removal 
times  and  independent  propabi 1 ities  p^,  i = 1,  2,  ...»  k,  for  removal  at  the 
it*1  inspection.  The  overall  emulative  distribution  function  of  engine 
removal  time  is  the  probability  that  removal  time  occurs  prior  to  or  at  age  t 

i (t) 

F(t)  = 1-  (1-F,(t))  n (1-P.)  (4.1) 

1 i=l  1 

where  i(t)  is  the  index  of  the  last  inspection  prior  to  or  at  age  t.  This 
engine  removal  time  model  assumes  that  inspections  occur  at  fixed  times  and 
that  inspection  removal  times  and  usage  removal  times  are  statistically  inde- 
pendent. The  latter  is  a testable  hypothesis  and  a likelihood  ratio  test  to 
do  so  has  been  prepared  but  not  programmed  or  tested.  The  former  assumption 
is  an  idealization  from  reality  that  may  be  an  adequate  representation  since 
engine  removal  codes  distinguish  inspection  removals.  If  not,  a similar  model 
could  be  constructed  to  represent  inspection  times  which  vary  around  a fixed 
value. 
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The  maximum  likelihood  estimator  of  F-|(t)  and  the  has  been  derived 
in  Appendix  A.  It  is  based  on  a sample  consisting  of  engine  usage  and 
inspection  removal  times  and  survivors'  ages.  When  combined  according  to 
formula  (4.1)  to  obtain  an  estimate  of  engine  removal  time  cumulative  distri- 
bution function  regardless  of  cause,  the  product  limit  estimator  of  Kaplan 
and  Meier  [4  ] is  obtained, 

F(t)  = l-n(N-r)/(N-r+l ) (4.2) 

r 

where  N is  the  total  sample  size  and  the  index  r ranges  over  the  indexes  of 
ordered  data  from  smallest  to  largest  including  in  the  product  above  only 
those  indexes  for  removal  times  prior  to  or  at  time  t.  An  example  of  the 
computation  of  this  estimator  is  in  Appendix  A. 

It  is  very  easy  to  convert  the  product  limit  estimator  to  an  estimator 
of  the  cumulative  failure  rate  function 

R(t)  = -log(l-F(t)  ).  (4.3) 


The  estimator  is  obtained  from  the  logarithm  of  (4.2) 


R(t)  = l (N-r)/(N-r+l ) 
r 

The  relation  between  actuarial  failure  rates  and  the  cumulative  failure 
rate  function  R(t)  is 

R(iAt)  = £ Ri 
j=l 

where  R^  is  the  true  interval  failure  rate  or  the  true  probability  of  engine 
removal  between  ages  (i-l)At  and  i At  given  survival  to  age  (i-l)At.  Conse- 
quently comparison  of  estimators  can  be  made  on  the  basis  of  interval  failure 
rates  or  cumulative  failure  rates,  or  even  on  cumulative  distribution  functions 
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since  the  actuarial  failure  rates  can  be  converted  to  a cumulative  distri- 
bution function  by 

i ( t ) - 1 

F ( t ) = 1-exp  (-1  R.  - Ri(t)  (t-(i-l)At)  ) 

j=l 

where  i(t)-l  is  the  index  of  the  actuarial  age  interval  containing  age  t. 

This  comparison  is  done  in  Section  5. 

Furthermore,  it  is  possible  to  use  each  estimator,  actuarial,  product 
limit  and  maximum  likelihood,  to  generate  engine  lives  and  construct  simulated 
operation  of  planes  for  comparison  of  how  well  the  estimators  predict  engine 
replacement  requirements.  This  comparison  is  made  in  Part  II,  Section  2. 
However,  it  is  a well  known  fact  that  functions  of  maximum  likelihood 
estimators  are  also  maximum  likelihood,  and  since  all  three  estimators  are 
maximum  likelihood  estimators  for  their  sample  information  and  statistical 
models  of  removal  times,  any  estimators  obtained  from  them  of  replacement 
requirements  are  maximum  likelihood  estimators.  Theoretical  comparison 
of  the  three  estimators  can  be  done  by  computing  their  information  matrices 
Wilks  [10]  Chapter  12. 

A computer  program  has  been  prepared  to  calculate  the  maximum  likelihood 
estimator  of  F^(t)  and  the  P.  i = 1 , 2,  . . . , k and  the  product  limit  estima- 
tors of  engine  removal  time,  cumulative  distribution  function,  and  failure 
rate  function.  (The  flow  chart  is  shown  in  Figure  2 and  the  procedure  is 
demonstrated  in  the  next  section). 

It  is  necessary  to  point  out  again  that  the  product  limit  estimators 
and  the  maximum  likelihood  estimators  are  based  on  sample  data  containing 
engine  removal  times  and  survivors'  ages,  not  the  information  on  when  the 
engines  began  operating  during  the  period  of  data  collection.  Further 
derivation  is  necessary  to  obtain  maximum  likelihood  estimators  based  on 
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sample  removal  times  and  ages  and  including  operating  time  experience  in 
a fixed  calendar  period.  It  is  also  worth  questioning  whether  it  would  be 
better  not  to  limit  experience  to  a fixed  calendar  period  for  some  purposes. 

It  seems  that  limitation  of  data  to  a particular  calendar  period  is  necessary 
only  when  deciding  whether  the  experience  during  that  period  differs  from 
earlier  experience.  If  the  decision  is  that  there  is  no  difference,  then 
the  product  limit  and  maximum  likelihood  estimators  may  be  used  to  obtain 
interval  failure  rates  or  any  of  the  other  reports  computed  from  actuarial 
failure  rates. 

One  last  remark  about  the  product  limit  and  maximum  likelihood  estima- 
tors computation  should  be  made.  They  both  require  some  sorting  of  data 
on  removal  times  and  survivors'  ages.  Modern  computers  can  sort  huge  amounts 
of  data  in  little  time,  but  the  maximum  likelihood  estimator  is  by  its 
nature  quicker  to  calculate  since  removal  times  and  survivors  ages  may  be 
sorted  separately.  The  product  limit  estimator  must  first  order  all  observa- 
tions. The  actuarial  estimator  requires  only  that  numbers  of  exposures  and 
deaths  per  age  interval  be  counted.  However,  the  flexible  actuarial  estimator 
proposed  in  Section  3 will  require  sorting. 


5.  Comparison  of  tne  Estimators 

, Estimation  of  simulated  engine  removal  ages  has  been  done  here  to  show 
what  the  estimators  do  but  in  a much  simpler  example  than  real  life. 

The  objective  is  to  simulate  removal  age  data,  compute  all  three  estimators 
of  the  removal  age  cumulative  distribution  function  and  failure  rates,  and 
compare  these  with  the  theoretical  cumulative  distribution  function 
used  to  simulate  the  removal  age  data. 

The  input  parameters  are  twenty  engines  with  starting  age  zero  and 
max  time,  two  hundred  hours.  There  will  be  one  inspection  at  one  hundred 
hours  and  the  probability  of  an  engine  is  removed  at  this  inspection  is  one- 
fifth.  The  probability  of  removal  at  the  maximum  time  of  200  hours  is  one 
to  insure  that  all  surviving  engines  are  removed  at  max  time. 

The  removal  age  data  was  generated  from  three  formulas.  The  usage 
removal  age  random  variable  is  generated  from  an  exponential  distribution 
with  a mean  removal  age  of  two  hundred  hours.  The  survivor's  age  is  a 
random  variable  Y generated  using  a uniform  distribution  from  age  zero  to 
four  hundred.  Four  hundred  hours  was  used  to  assure  that  there  would  not 
be  too  many  survivors  up  to  the  maximum  operating  time  of  two  hundred 
hours.  The  inspection  removal  time  random  variable  Z is  generated  for  one 
hundred  and  two  hundred  hours  by  the  probabil ities  given  earlier. 

The  engine  removal  age  or  survivors'  age  is  the  minimum  of  these  three 
variables. 

The  resulting  random  variable,min  (X,Y,Z),  has  the  theoretical  cumulative 
distribution  function 
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T 


i 


i(t) 

1 . e‘t/20°  • (l-t/400)  ^ (1-Pi)  (5.1) 


where 


i(t) 


0 t < 100  hours  prior  to  inspection 

1 1 00<_  t < 200  hours 

2 t > 200  hours. 


(If  i(t)  = 0,  the  product  is  1.0 

The  objective  is  to  estimate  the  engine  removal  age  min  (X,Z)  which 
is  the  minimum  of  usage  removal  age  or  inspection  removal  time  (including 
max  time  as  the  second  inspection)  in  the  presence  of  the  nuisance  variable 
Y,  survivor's  age. 

The  theoretical  cumulative  distribution  function  of  min  (X,Z)  the 
engine  removal  age  is 

0 t < 0 


F(t)  = i-e_t/20° 


i (t) 

n (1-Pi)  = 
i = l 


ll  - e’t/20° 
|1  - .Se'^200 
1 


0 < t < 100 
100  i t < 200 
t > 200 


(5.2) 


which  is  plotted  in  Figure  3.  The  maximum  likelihood  estimator  also  gives 
separate  estimates  of  the  cumulative  distribution  function  of  usage  removal 
age  and  the  inspection  removal  probabilities  for  use  in  engine  diagnosis 
and  prediction  of  replacement  requirements  due  to  inspection  removals. 

The  construction  of  Table  1 is  as  follows:  Column  1 is  the  index 

of  engines  available,  for  this  example  that  is  twenty  engines.  These  could 
be  the  engine  serial  numbers  for  real  data.  Columns  2,  4,  and  6 are 
random  numbers  taken  from  Mihram  [ 2 ].  A different  random  number  is 
generated  for  each  engine's  age  at  usage  removal,  survivor's  age,  and 
inspection  removal  times.  Column  3 uses  the  inverse  transformation  method 
(Fishman  [ 1 ])  for  generating  exponentially  distribution  random  variables. 
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TABLE  1 

SIMULATED  ENGINE  REMOVAL  AGF  DATA 


Col umn 
Number 

1 

2 

3 

4 

5 

hours 

6 

7 

8 

Column 

Number 

Engine 

Number 

Random 

aNumbers 

X=-L0G(RN)*200 

RN's 

Y=RN*400 

RN's 

Z 

b 

MIN(X,Y,Z) 

1 

.71264 

68 

.362 

145 

.295 

200 

68 

2 

.84665 

32 

.469 

188 

.555 

200 

32 

3 

.24364 

282 

.015 

6 

.341 

200 

e> 

4 

.96418 

7 

.902 

361 

.307 

200 

7 

5 

.21437 

308 

.541 

216 

.634 

200 

200 

6 

.95026 

10 

.957 

383 

.088 

100 

10 

7 

.02731 

720 

.066 

26 

.519 

100 

@ 

8 

.24453 

282 

.631 

252 

.474 

200 

200 

9 

. 1 5006 

379 

.812 

325 

.933 

200 

200 

10 

.35218 

209 

.581 

232 

.077 

100 

100 

11 

.44093 

164 

.497 

199 

.513 

200 

164 

12 

.54629 

121 

.25034 

100.1 

.788 

200 

100 

13 

.66392 

82 

.729 

292 

.907 

200 

82 

14 

.04394 

625 

.326 

130 

.586 

2 00 

© 

15 

.28697 

250 

.444 

178 

.386 

200 

© 

16 

.95000 

10 

.263 

105 

.451 

200 

10 

17 

.15711 

370 

.246 

98 

.175 

100 

18 

.56211 

115 

.102 

41 

.708 

200 

© 

19 

.14635 

384 

.568 

227 

.477 

200 

200 

20 

.63087 

92 

.713 

285 

.489 

200 

92 

a.  Mihram  [ 2 ] Random  Number  Tables  p.  499. 

b.  Survivor's  ages  are  circled. 


Column  5 uses  a uniform  distribution  and  rescales  the  random  number  to  the 
interval  0-400  hours.  Column  7 generates  inspection  or  max  time  removals. 

If  the  random  number  is  less  than  0.2,  Z is  an  insoection  time  of  one 
hundred  hours,  otherwise  Z is  the  max  time  of  two  hundred  hours. 

Column  8 is  the  minimum  time  of  columns  3,  5 or  7.  These  ages  constitute 
the  data  of  a simulation  of  twenty  engines  operating  to  a possible  max  time  of 
two  hundred  hours. 

The  removal  ages  and  survivors  ages  from  Column  8 are  used  to  compute 
the  actuarial  failure  rates  in  Table  2.  In  this  table,  there  are  five 
actuarial  age  intervals  assumed  to  have  width  forty  hours,  and  all  maximum 
time  removals  are  collected  as  a sixth  interval  of  zero  width.  Column  1 is  the 
interval  index.  Column  2 gives  the  age  span  for  each  of  the  intervals. 

Column  3 is  the  number  of  removals  during  each  interval.  These  are 
obtained  from  the  removal  times  in  Table  1.,  Col.  8.  For  instance,  engines 
number  2,  4,  6,  and  16  are  removed  at  times  33,  7,  10  and  10  respectively 
so  there  are  four  removals  in  interval  one.  Column  4 is  the  number  of 
survivors'  ages  contained  in  each  interval.  These  are  found  in  the  same 
way  as  the  removals  times  in  column  3.  Column  5 gives  the  number  of  full 
exposures  for  each  interval.  This  is  the  number  of  engines  that  survived 
to  an  age  greater  than  the  end  of  the  interval  plus  the  number  removed  during 
the  interval.  The  number  of  partial  exposures  is  the  proportion  of  an  age 
interval  during  which  engines  recorded  as  survivors'  ages  operated.  For 
instance,  during  interval  one,  engine  3 operated  6 hours  or  .15  of  the 
interval  and  engine  7 operated  .65  of  the  interval  for  a total  of  .8  partial 
exposure^  in  interval  one. 
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TABLE  2 


i 


COMPUTATION 

OF  ACTUARIAL 

FAILURE 

RATES 

Column 

Number 

Column 

Name 

1 

Interval 

2 

Span 

3 

^Removals 

4 

#Survi  vor 

5 

Full 

^Exposures 

6 

Partial 

Exposures 

7 

^Exposures 

8 

Fa i 1 ure 
Rate 

1 

0-40 

4 

2 

18 

.8 

18.8 

.213 

1 

2 

40-80 

1 

1 

13 

+ 

.025 

13.025 

.077 

3 

80-120 

4 

1 

11 

+ 

.45 

11.45 

.349 

4 

120-160 

0 

1 

6 

+ 

.25 

6.25 

0 

5 

160-200 

1 

1 

5 

+ 

.45 

5.45 

.32 

max  time 

200 

_4 

_ 

4 

4 

1 

(i 

14 

6 

20  engines  total 

TABLE  3 

COMPARISON  OF  ACTUARIAL  ESTIMATOR  AND  THEORETICAL  DISTRIBUTION 


lol  umn 
Number 
lol  umn 
Name 

1 

Interval 

2 

Ri 

Failure 

Rate 

3 

cdf.  end 
of  Interval 

4 

Survi vor 
Prob. 
Il-F(t)] 

5 

Theoretical 
c.d.f.  of 
Removal  Times 

6 

Theoretical 
Interval 
Failure  Rate 

1 

(40) 

.213 

.192 

.808 

.181 

- 

.2 

2 

(80) 

.077 

.252 

.748 

.330 

- 

.2 

3 

(120) 

.349 

.472 

.528 

.561 

- 

.483 

4 

(160) 

0 

.472 

.528 

.641 

- 

.2 

5 

(200') 

.184 

.561 

.439 

.706 

- 

.2 

max 

time 

MT 

(200) 

1 

1 

0 

1 .0 

1 .0 

I 
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Column  7 is  the  total  number  of  exposures  which  is  the  sum  of  Columns  5 and  6. 
Column  8 is  the  failure  rate  for  each  interval  which  is  obtained  by  dividing 
the  number  of  removals  (Column  3)  by  the  number  of  exposures  (Column  7). 

In  Table  3,  there  is  a comparison  of  the  cumulative  distribution  function 
and  failure  rates  for  the  actuarial  estimator  and  theoretical  distribution 
(5.2).  Column  1 gives  the  number  of  the  interval.  Column  2 gives  the 
actuarial  failure  rate,  R ^ which  is  column  8 of  Table  2.  Column  3 is  the 
actuarial  cumulative  distribution  function,  F(t^),  where  t^  is  the  end  of 
the  interval.  This  function  is  obtained  from 

F.  (t-)  1 - exp  [-  I R.] 

A j = l J 

where  R.  is  the  failure  rate  in  column  2.  The  survivor  probabilities  in 

s) 

column  4 are  the  complements  of  column  3,  l-F(t-j).  Column  5 gives  the 
theoretical  cumulative  distribution  function  of  removal  times  by  the 
equation  (5.2)  where  t is  evaluated  at  the  end  of  each  interval.  To  find 
the  theoretical  interval  failure  rates,  column  5 is  used  and  the  theoretical 
interval  failure  rate  function  R(t.j ) is  found  solving  iteratively 

i 

F(t.)  = 1 - exp  [-  £ R ( t • ) ] i = 1,2,. ..,6 

1 j=l  J 

for  R(t^)  in  each  interval  where  t^  = 40,  80,  120,  180,  200  and  200. 

Figure  2 shows  the  comparison  of  the  actuarial  and  theoretical  interval 
failure  rates  as  a graphical  representation.  Figure  3 shows  the  comparison 
of  the  actuarial  and  theoretical  cumulative  distribution  functions. 

It  should  be  recognized  that  deviation  of  the  actuarial  estimator  from  the 
theoretical  value  is  due  the  very  small  sample  size  of  14  observed  removal 
times  and  the  fact  that  the  data  was  simulated  from  the  theoretical 
distribution  function  (5.1). 
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The  product  limit  estimator  of  the  cumulative  distribution  function 
and  failure  rates  is  calculated  from  the  data  in  Table  1 as  follows: 

Columns  1 and  3 are  the  same  as  columns  1 and  8 in  Table  1 except  that  the 
removal  times  are  reordered  along  with  their  indexes.  Column  2 is  the  re- 
moval code  useful  here  and  in  calculation  of  the  maximum  likelihood 
estimator.  Column  4 is  the  product  limit  estimate  of  the  percent  surviving 
or  tail  cumulative  distribution  function  at  the  time  or  age  values  in 
column  3.  The  formula  as  given  in  Kaplan  and  Meier  [ 4 J is 

II  (N-r)/(N-r+l ) 

r 

where  N is  the  sample  size  20  and  the  index  r ranges  from  1 to  20  taking 
on  only  the  values  for  inspection  and  usage  removals  in  column  3.  Column  5 
is  the  complement  of  column  4,  the  product  limit  estimator  of  the  cumulative 
distribution  function  (5.2).  Column  6 is  the  value  of  the  failure  rate 
function  assuming  that  it  is  constant  during  the  interval.  It  is  the  log 
of  (N-r)/(N-r+l ) for  intervals  between  engine  removal  age  r and  the  next 
subsequent  engine  removal.  Column  7 is  the  cumulative  failure  rate  function 
at  the  age  in  column  3,  which  is  also  the  natural  logarithm  of  Column  4. 

The  product  limit  estimator  and  the  theoretical  cumulative  distribu- 
tion functions  are  plotted  in  figure  4 and  their  corresponding  cumulative 
failure  rates  are  in  figure  5.  For  comparison  of  the  product  limit  estimator 
with  the  actuarial  estimator  of  failure  rates,  the  product  limit  estimator 
of  failure  rates  during  each  40  hour  age  interval  must  be  calculated  from 
the  cumulative  failure  rate  estimator.  The  comparison  of  product  limit 
estimator  with  the  theoretical  interval  failure  rates  is  shown  in  figure  6. 
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TA3LE  4 


Column 
Number 
Col umn 
Name 


PRODUCT  LIMIT  CU'1ULATIVE  FAILURE  RATE  FUNCTION 
1 2 3 4 5 6 7 


Index 

Code 

Age 

ll(N-r)/(N-r+l ) 
r 

FPL ' fc) 

ri* 

R(t) 

1 

Survi  vor 

0 

- 

.0 

2 

Usage 

7 

.947 

.053 

.054 

.054 

3 

Usage 

10 

.895 

.105 

.057 

.111 

4 

Usage 

10 

.842 

.158 

.061 

.172 

5 

Survivor 

© 

- 

.061 

- 

6 

Usage 

32 

.786 

.214 

.069 

.241 

7 

Survi vor 

© 

- 

.064 

- 

8 

Usage 

68 

.726 

.274 

.080 

.321 

9 

Usage 

82 

.665 

.335 

.087 

.408 

10 

Usage 

92 

.605 

.395 

.095 

.503 

11 

Survivor 

© 

- 

.095 

- 

12 

Inspection 

100 

.537 

.453 

.118 

.621 

13 

Inspection 

100 

.470 

.530 

.134 

.755 

14 

Survivor 

© 

- 

.134 

- 

15 

Usage 

164 

.392 

.608 

.182 

.927 

16 

Survi vor 

© 

- 

.182 

- 

17 

Max  time 

200 

.294 

.706 

.288 

1.225 

18 

Max  time 

200 

.405 

1.630 

19 

Max  time 

200 

.693 

2.32 

20 

Max  time 

200 

0 

1.0 

oo 

00 

*Failure  rate  function  is  assumed  constant  between  observations  of  actual 

removals . 


.82 

.93 

1.73 
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In  comparison  with  actuarial  failure  rates  in  Figure  2, 
the  product  limit  estimator  is  closer  than  the  actuarial  estimator  to  the 
theoretical  interval  failure  rate  except  in  intervals  1 and  5. 
expected  that  this  superiority  will  increase  as  sample  size  grows 
because  the  product  limit  estimator  uses  the  actual  values  of  removal  times 
and  survivors'  ages  instead  of  lumping  the  time  values  into  interval  counts 
exposures  and  removals. 

It  is  very  easy  to  calculate  interval  failure  rates  for  any  width 
of  age  intervals  from  the  cumulative  failure  rate  function  in  column  8 
or  any  cumulative  failure  rate  function  R(t),  t ^ 0.  If  At  is  the  desired 
width  of  age  interval,  then  the  interval  failure  rate  here  denoted  as 


Rpli = R(i At)  - R ( ( i - i ) At) . 

For  example,  in  order  to  obtain  the  product  limit  estimator  of  the  failure 
rate  in  interval  2,  40  to  80  hours,  use  linear  interpolation  between 
usage  removals  at  33  and  68  hours  to  obtain  R(40)  = .25  and  between  usage 
removals  at  68  and  82  hours  to  obtain  R(80)  = .38.  (These  values  were  read 
from  figure  5).  The  interval  failure  rate  Rp^  is  the  difference,  .13. 

The  maximum  likelihood  estimators  for  the  usage  removal  time  cumulative 
distribution  function  and  the  probability  of  removal  at  inspection  have 
also  been  computed  according  to  formulas  in  Appendix  A.  (The  maximum 
likelihood  estimator  of  the  removal  time  cumulative  distribution  function 
regardless  of  whether  it  is  an  inspection  or  usage  removal  is  the  product 
limit  estimator.)  The  estimate  of  the  probability  of  removal  at  the  first 
inspection  is 

number  removed  at  inspection  _ 2 _ ^2 

number  surviving  to  inspection  9 
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(the  theoretical  value  was  0.2)  and  the  estimate  of  the  usage  removal  time 
cumulative  distribution  function  is  equivalent  to  Figure  4 without  the 
jump  at  100  hours.  This  is  plotted  against  the  theoretical  formula  used 
to  generate  usage  removal  times  in  Figure  7.  The  jump  at  200  hours 
represents  the  probability  of  surviving  usage  removal  to  max  time. 

A computer  program  for  comparison  of  estimators  on  several  bases  has 
been  written.  In  the  remainder  of  this  section,  estimators  are  compared  on 
the  basis  of  their  cumulative  distribution  function.  In  Part  II,  Section  8, 
estimators  are  compared  on  the  basis  of  how  well  they  predict  replacement 
requirements.  The  flow  chart  for  computation  of  the  maximum  likelihood 
estimators  is  shown  in  Figure  8.  In  Figure  9 the  actuarial  and  product 
limit  estimators  are  compared  with  the  theoretical  distribution  with  three 
inspections  at  50,  100  and  150  hours.  No  survivors'  ages  were  generated 
(no  censoring),  but  the  rest  of  the  program  generated  data  and  computed 
estimators  as  was  done  earlier  in  this  section.  Figures  10  and  11  are  the 
maximum  likelihood  estimators  of  the  usage  removal  time  cumulative  dis- 
tribution function  which  was  calculated  from  data  generated  from  an 
exponential  usage  removal  time  with  mean  133.5  and  inspection  at  100  hours. 
Figure  11  shows  the  estimator  in  the  presence  of  some  survivors'  ages. 

Figure  12  shows  the  cumulative  failure  rate  function  from  figure  10 
plotted  along  with  the  theoretical  failure  rate  function.  The  slight  low 
side  bias  in  all  three  figures  is  due  to  the  fact  that  the  simulated  data 
contained  fewer  usage  removals  than  expected,  not  because  of  systematic 
bias  in  the  estimators.  A best  fit  linear  regression  to  the  cumulative 

failure  rate  estimator  plotted  in  Figure  12  gave  a slope  of  .0070  and  an 

’ ? 

intercept  of  -.0064  with  an  R value  of  .9966  indicating  the  data  almost 

| 33 

d 
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certainly  come  from  a linear  model  for  a cumulative  failure  rate  function 
which  it  did  since  the  cumulative  failure  rate  function  for  exponential 
usage  removal  ages  is  linear. 
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Figure  8,  Maximum 
Likelihood  Estimator 
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6.  Updating  the  Data  Set  for  Computing  Official  Failure  Rat^ 

There  are  several  tests  used  to  infer  whether  engine  performance 
is  changing.  (All  of  these  are  found  in  AFLCM  66-17  [11]).  The  "Control 
Factor"  test  uses  the  normed  difference  between  actual  and  expected  failures 
experienced  in  an  actuarial  age  interval.  If  the  magnitude  of  the  control 
factor  is  too  large,  the  official  failure  rate  used  to  estimate  expected 
failures  is  suspect  (or  perhaps  usage  by  a base  or  squadron  differs 
significantly  from  the  official  failure  rate).  This  is  method  A section  9 
of  AFLCM  66-17.  The  "statistical  test"  (page  4-10)  is  computed  to  compare 
crude  and  official  failure  rates  (or  crude  and  smoothed  failure  rates)  on 
an  interval  by  interval  basis.  If  the  interval  failure  rate  is  significantly 
different  from  the  official  failure  rate,  it  may  be  changed  proportional  to 
the  ratio  of  the  official  actuarial  life  expectancy  and  the  actual  life 
expectancy.  The  third  test,  "Method  D"  of  section  9,  monitors  flying  hours 
per  failure  over  several  calendar  time  periods  with  a control  chart. 

Only  in  method  D is  inference  based  on  the  experience  of  an  entire 
population  of  one  type  of  engines.  The  other  two  tests  base  inference  on 
test  statistics  calculated  for  each  actuarial  age  interval,  and  consequently 
the  conclusion  for  adjacent  age  intervals  may  be  contradictory.  These  two 
tests  may  be  very  sensitive  for  engine  diagnosis  to  changes  in  engine 
component  performance,  but  they  may  not  be  the  appropriate  way  to  determine 
whether  replacement  requirements  will  change. 

The  general  questions  that  seem  to  be  addressed  in  the  statistical 
tests  are  whether  the  estimates  of  engine  performance  have  changed  and  whether 
engines  seem  to  be  failing  at  a particularly  high  rate  during  some  age 
interval.  There  are  two  purposes  here,  engine  diagnosis  and  replacement 


*1 

requirements  prediction.  Engine  diagnosis  can  be  assisted  by  determining 
whether  there  is  a significant  number  of  failures  near  some  age.  If  so 
then  there  may  be  an  engine  component  whose  operating  life  is  significantly 
less  than  the  engine  life  itself  and  the  component  should  be  replaced  prior  to  its 
failure.  The  two  tests,  "control  factor"  and  "statistical  test"  have  their 
discriminatory  power  reduced  by  the  smoothing  technique  mentioned  in  section 
2 of  this  report.  Smoothing  tends  to  smooth  out  peaks  in  failure  rates  • 
that  may  be  caused  by  failure  of  some  engine  component.  Those  two  sta- 
tistical tests  will  have  to  be  revaluated  if  smoothing  is  stopped  and 
more  sensitive  engine  diagnosis  is  possible  using  the  estimators  proposed 
in  section  4 . 

The  remainder  of  this  section  is  to  determine  whether  the  official 
failure  rates  are  adequate  for  prediction  of  replacement  requirements. 

Future  engine  requirements  depend  on  the  entire  set  of  failure  rates  through- 
out the  operating  time  from  new  to  max  time,  so  comparison  of  the  official 
failure  rates  to  current  failure  rates  should  be  made  on  an  aggregate  basis. 

It  is  advocated  that  the  comparison  should  be  made  of  the  original  data 
instead  of  the  failure  rates  computed  from  the  data  for  the  same  reason 
in  section  1 , that  the  actuarial  failure  rates  are  based  on  a summary 
of  the  data.  The  engine  life  data  itself  should  be  used  to  answer  questions 
about  engine  lives  and  consequent  replacement  requirements. 

There  are  several  ways  to  statistically  formulate  the  question  of 
whether  the  official  failure  rates  still  indicate  engine  performance. 

1.  Does  the  current  calendar  period  of  engine  experience  resemble 
the  experience  used  to  compute  official  failure  rates?  Here  experience 
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refers  to  recorded  engine  removal  times,  and  exposures  during  the 
calendar  periods  being  compared. 

2.  Do  the  recorded  ages  at  removal  during  the  current  calendar 
period  and  the  ages  of  engines  still  operating  at  the  end  of  this 
period  come  from  the  same  population  of  engine  lives  as  the  corresponding 
data  used  to  calculate  the  official  failure  rates? 

These  questions  differ  because  an  engine  removed  during  the  current 
calendar  period  may  have  acquired  a significant  proportion  of  its  operation 
in  prior  calendar  periods.  A statistical  procedure  for  answering  question  2 
is  provided  in  this  section, and  a method  for  modifying  this  procedure  to 
answer  question  1 will  be  described.  The  statistical  procedure  computes 
a measure  of  how  likely  it  is  that  the  current  calendar  period's  data  could 
have  come  from  the  same  population  as  past  data.  This  measure  is  computed 
for  every  past  calendar  period  for  which  data  has  been  collected  taking 
continguous  periods  ending  with  the  beginning  of  the  current  period. 

The  Aerospace  Engine  Life  Committee  may  then  use  these  measures  and  their 
knowledge  of  the  actual  operating  conditions  during  the  current  period 
to  decide  whether  to  update  the  data  set  to  include  current  data  and 
whether  to  eliminate  the  oldest  data. 

The  test  uses  a Wilcoxon  type  rank  test  suggested  by  Gehan  [6]. 

The  Gehan  test  statistic  was  used  without  the  later  modification  by 
Efron  [19].  A flowchart  of  this  test  is  in  Figure  13.  The  performance 
of  this  test  appears  to  be  as  good  as  any  two  sample  test  with  the  possible 
exception  of  the  Kolmogorov- Smirnov  type  test.  (Several  two  sample  tests 
were  compared  in  Lee,  Desu  and  Gehan  [17]  and  Gehan  and  Thomas  [ lb].) 

The  reason  for  the  exception  of  the  Kolmogorov-Smirnov  type  test  is  that 
it  has  not  yet  been  adapted  to  incorporate  censoring  such  as  experienced 
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when  sample  information  includes  the  ages  of  surviving  engines.  This  can 
be  done  by  extension  of  the  truncated  Kolmogorov-Smirnov  two  sample  test 
in  Barr  and  Davidson  [20]  and  Koziol  and  Byar  [21  ]. 

The  Wilcoxon  type  rank  test  used  in  the  updating  program  adds  a 
plus  or  minus  one  to  the  test  statistic  W if  it  can  be  concluded  that  a 
engine  age  at  removal  from  the  current  period  is  or  will  be  larger  than  an 
engine  age  at  removal  in  the  past.  The  proposed  modification  to  this  test 
to  answer  question  1,  whether  the  current  period's  experience  is  the  same 
as  the  past,  is  to  weight  the  plus  or  minus  one  by  the  proportion  of  an 
engines  exposure  that  is  in  the  current  period.  This  modification  has  not 
been  evaluated. 
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Figure  13 

Update  Main  Program 


SECTION  II 

PREDICTION  OF  REPLACEMENT  REQUIREMENTS 
7.  Prediction  of  Engine  Requirements  from  the  Actuarial  Failure  Rates 

Two  fundamental  ingredients  are  necessary  to  predict  replacement 
requi rements,  engine  removals  which  require  an  engine  from  spares  stock  to 
replace  the  one  that  was  removed;  some  probability  measure  of  the  amount 
of  operating  till  removal  must  be  mixed  with  the  amount  of  operating 
time  demanded  from  each  aircraft  to  obtain  measures  of  replacement 
requirements.  There  are  several  measures  of  replacement  requirements  that 
are  needed  for  engine  management: 

1.  The  expected  numbers  of  engine  repairs  and  rebuilds 

are  needed  for  spare  parts  provisioning  and  work  load  planning. 

2.  Upper  confidence  intervals  or  quantiles  of  engine 
requirements  are  needed  to  determine  the  number  of  spare 
engines  necessary  to  meet  mission  requirements  with  a 
sufficiently  high  probability. 

3.  Predictions  of  the  actual  time  when  operating  engines 
will  need  inspection,  repair,  or  rebuilding  are  necessary 
for  workload  and  maintenance  planning. 

The  actuarial  method  contains  elements  designed  to  meet  each  of  these 
requi rements . 

The  "actuarial  engine  life"  is  an  estimate  of  the  time  till  removal  of 
a new  or  rebuilt  engine.  If  the  overhaul  failure  rates  are  used  in  the 
computation, then  the  actuarial  engine  life  is  an  estimate  of  the  expected 
time  till  overhaul.  If  base  maintenance  removal  rates  are  used, then  the 
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actuarial  engine  life  is  an  estimate  of  the  time  till  first  removal  for 
base  repair  starting  with  a new  or  rebuilt  engine.  The  expected  time 
till  removal  could  be  converted  into  an  estimate  of  demand  requirements 
except  that  not  every  engine  is  new  at  the  beginning  of  the  period  for 
which  replacement  predictions  are  needed.  The  actuarial  method  includes 
a procedure  for  "operating"  on  paper  the  entire  inventory  through  a calendar 
period  to  determine  replacement  requirements  from  the  input  information  on 
inventory,  ages,  flying  hour  program  and  actuarial  failure  rates.  This 
procedure  is  described  in  T.O.  00-25-128  Chapter  IX  [7],  T.O.  00-25-217 
Chapter  IX  [8]  and  AFLC  Manual  66-17  Chapter  6 [11].  It  is  equivalent  to 
computing  the  conditional  expected  operating  time  till  replacement  given 
engine  age  at  the  beginning  of  the  calendar  period  under  consideration. 

This  expected  residual  operating  time  is  in  turn  used  to  estimate  replace- 
ment requirements  to  meet  the  flying  hour  program  by  assuming  the  flying 
hours  are  equally  allocated  among  all  aircraft  and  each  engine  contributes 
its  expected  residual  life  to  meet  the  program. 

This  procedure  for  predicting  expected  replacement  requirements  is 
based  on  two  fallacies.  For  example,  if  every  engine  is  new  and  the 
expected  operating  time  till  replacement  is  300  hours,  then  it  is  natural 
and  convenient  to  conlcude  that  to  fly  1000  hours  will  require  two  replace- 
ments or  three  engines  total.  In  general  this  procedure  is  as  follows:  let 

ET  denote  the  expected  time  till  replacement,  t the  flying  hours  required 
from  a single  engine  plane,  and  N(t)  the  number  of  engines  required  to  meet 
the  flying  hour  requirement  of  t hours.  An  obvious  estimate  of  N(t)  is 
N(t)  = t/ET. 

Unfortunately  this  estimate  is  unbiased  only  in  the  very  special  case  of 
a single  engine  plane  with  new  engines  and  only  when  the  engine  lives  are 
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independent  identically  distributed  exponential  random  variables  (Cox 
Chapter  2 [12]).  This  is  fallacy  #1.  It  is  faintly  possible  that  the 
procedure  for  estimating  replacement  requirements  is  asymptotically  un- 
biased for  a large  fleet  of  simultaneously  operating  aircraft.  However, 
one  condition  for  this  to  be  true  is  that  all  engines  have  identically  dis- 
tributed operating  times  until  reolacement  (Cox  Chapter  6 [12]).  This 
assumption  is  not  true  for  engines  that  have  accumulated  some  operating 
hours  at  the  beginning  of  the  calendar  period  under  consideration  unless 
all  engines  have  exponentially  distributed  times  till  replacement.  Some 
indication  of  the  error  and  direction  of  the  bias  in  the  estimate  of  N(t) 
is  possible  from  analytical  results.  This  is  discussed  further  in  section  9, 
part  II.  Meanwhile,  the  only  unbiased  estimator  of  the  expected  engine  re- 
quirements must  be  obtained  by  simulation. 

A procedure  for  computing  the  number  of  spare  engines  needed  to  meet 
flying  hour  requirements  with  90 % probability  is  given  in  AFM  400-1  Chapter 
8 [13].  The  procedure  is  based  on  the  approximation  alluded  to  in  the  pre- 
vious paragraph.  The  approximation  is  that  replacement  requirements  will 
have  a Poisson  distribution  with  a parameter  which  is  the  inverse  of  the 
mean  time  between  replacements.  The  90th  Dercentile  of  the  Poisson  distribu- 
tion is  used  to  determine  replacement  requirements.  Again,  this  approximation 
is  not  valid  because  engine  ages  at  replacement  are  not  exponential,  but  it 
may  be  possible  to  strengthen  the  approximation  or  change  to  an  alternative 
normal  approximation.  This  is  considered  in  Section  3,  Part  II.  Meanwhile, 
simulation  remains  the  only  unbiased  procedure  for  predicting  confidence 
limits  on  spare  engine  requirements. 


The  second  fallacy  is  to  use  the  age  of  an  engine  at  the  beginning  of 
the  current  calendar  quarter  to  compute  the  residual  operating  time  to  removal 
as  if  the  engine  had  survived  to  the  age  recorded  at  the  beginning  of  the 
quarter.  Some  engines  used  during  a quarter  had  previously  failed,  and  their 
ages  at  the  beginning  of  the  quarter  may  be  their  age  at  failure  rather  than 
the  survival  time.  The  question  of  whether  this  fallacy  is  biasing  predictions 
of  engine  operating  times  requires  further  testing.  Meanwhile,  we  will  use 
the  engine  age  at  the  beginning  of  a quarter  as  a survivor's  age  rather  than 
a failure  time  to  keep  results  comparable  with  Air  Force  results. 

Prediction  of  the  third  measure  needed  for  engine  management,  the 
calendar  times  of  replacement  of  particular  engines,  has  not  been  done  by 
the  actuarial  method  although  the  ingredients  are  available  to  do  that.  It 
is  possible  to  do  this  by  simulation. 

L.A.  Coco  and  Dale  Plank  of  AFLC  have  written  computer  programs  to 
simulate  engine  requirements  based  on  the  flying  hour  program,  the  actuarial 
failure  rates,  and  the  ages  of  available  engines.  A simplified  flow 
chart  of  this  program  is  shown  in  Figure  14.  It  is  essentially  a randomized 
version  of  the  actuarial  method  for  estimating  replacement  requirements. 

Instead  of  estimating  the  proportion  of  removals  for  engines  in  each 
actuarial  age  interval  as  the  actuarial  failure  rate,  the  program 

simulates  whether  each  engine  survives  each  actuarial  age  interval  required 

of  it  to  achieve  the  flying  hour  program.  The  simulation  is  done  by 
generating  a random  decimal  between  zero  and  1.0.  If  the  decimal  is 

larger  than  the  failure  rate,  the  engine  survives  the  age  interval  and 

contributes  that  amount  of  time  to  the  flying  hour  program.  Since 
failure  rates  are  small,  many  random  decimals  are  needed  and  the  program 
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execution  time  is  so  long  that  only  one  run  through  the  flying  hour  program 
is  possible.  This  procedure  yields  an  unbiased  estimate  of  replacement 
requirements  but  no  indication  of  its  accuracy  and  no  way  to  predict  con- 
fidence limits  on  spare  engine  requirements.  Repeated  runs  through  the 
program  would  generate  many  values  of  replacement  requirements  for  each 
calendar  period  in  the  flying  hour  program  from  which  an  indication  of 
the  accuracy  of  estimated  replacements  and  of  the  confidence  limits  could 
be  obtained.  However,  the  alternative  simulation  approach  in  the  next 
section  is  more  efficient  and  is  not  dependent  on  the  type  of  estimator 
used  to  simulate  engine  age  at  replacement. 


I /NPvT  FA/lUAF  , 
/ HA7SS  AA/£>  I 
IFLY/a/O  Hoc/ A.  I 

FA  ooAA*r  I 


IrYT/AL  tZF  AQC  /<VT«Wi 
Co  OH  re  A X si  f f^o/.vf 

COVAJ  TEA  Tro 


ywtERMiNTK 
/wHtlHEREAfftME  TV 
-S^Rv  Mes  Ace  /hteamac  1^ 
IS  R^MDOM  HUM&eO.  cuss 

THAN  a CT<S4RJAt_  >/ 

F/ntuum  r fats' 

\ H / 


rfpiacc  o*j£ 
ptoxe  f/vo/me 

AMP  COUH  T rMA 
iast  /msPvAcs 
OPFMTtOH 
TtfuHAD  THF 
FLYIhG-  HcoA 
RESkOtKF  HfCAJlS 


C OUMF  THS  i AST 
/A/TfFVAL  TourAAb 
THC  FCYM&tiOO* 
REaotAf  Htenrs 


^fOtuATMAl 

Nt/Alf/ 


It  i 

Tc  Jfl 


X - I -t-1 


Figure  14.  Simulation  of 
Single  Engine  Aircraft 
Engine  Requirements. 
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a.  Simulation  of  Engine  Requirements  by  the  Next  Event  Method 

There  are  two  fundamentally  different  methods  for  simulating  time 
dependent  random  processes,  the  time  slice  method  and  the  next  event  method. 
The  simulation  program  described  in  the  previous  section  is  an  example  of 
the  time  slice  method  where  time  (in  that  program  time  is  engine  operating 
time)  is  advanced  in  fixed  increments  equal  to  the  width  of  an  actuarial  age 
interval.  The  alternative  method  advances  the  simulation  time  clock  to  the 
time  of  the  next  event,  which  for  the  engine  replacement  process  would  be  a 
replacement  or  the  receipt  cf  an  engine  from  base  repair  or  depot  to  be 
returned  to  stock.  In  both  simulation  methods,  it  is  possible  to  synchronize 
operating  time  and  calendar  time  so  that  returns  to  stock  and  demands 
for  replacements  may  be  expressed  in  calendar  time  or  in  operating  time. 

There  is  no  inherent  superiority  of  one  method  over  the  other.  The  time 
slice  method  requires  many  random  number  generations, but  the  next  event 
method  requires  a more  complex  time  advancement  procedure.  In  engine 
replacement  simulation,  the  time  slice  method  takes  so  much  computer  time 
that  repeat  runs  are  not  worthwhile,  and  the  output,  replacement  require- 
ments, is  an  estimate  with  no  indication  of  its  accuracy.  The  next  event 
method  quickly  repeats  many  simulations  of  a flying  hour  program  providing 
much  more  information  about  replacement  requirements. 

The  simulation  program  submitted  with  this  report  differs  in  several 
ways  from  the  programs  written  by  L.A.  Coco  and  Dale  Plank.  The  program 
here  simulates  one  random  removal  age  and  randomly  selects  whether  the 
removed  engine  receives  base  repair  or  rebuild  according  to  the  return 
ratio.  It  could  be  modified  to  generate  both  a repair  removal  age  and  a 
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rebuild  removal  age  and  direct  the  engine  to  base  repair  or  depot  depending 
whether  the  repair  removal  age  or  rebuild  removal  age  was  smaller.  On  the 
other  hand,  a replacement  model  will  be  described  in  section  9,  Part  II  which 
generates  combined  removals  and  splits  removed  engines  for  repair  and 
rebuild  according  to  an  age  dependent  probability  of  rebuild.  This  model 
may  be  appropriate  since  older  engines  are  more  likely  to  be  rebuilt.  The 
age  dependent  rebuild  probability  could  be  incorporated  into  the  computer 
program.  It  is  not  clear  at  this  time  which  option  is  preferable. 

Another  difference  of  the  program  submitted  with  this  report  is  that 
no  attempt  has  been  made  to  accurately  represent  the  repair  and  rebuild  time. 
Arbitrary  calendar  time  intervals  were  chosen  after  which  an  engine  removed 
for  repair  or  rebuild  rejoins  the  spare  engine  inventory.  Its  age  is  reset 
to  zero  if  it  was  rebuilt.  These  repair  and  rebuild  times  could  be 
randomized  if  desired,  but  AFLC  already  has  a computer  program,  the  JEMS 
(Jet  Engine  Management  Simulator),  which  thoroughly  represents  every 
aspect  of  repair  and  rebuilding.  However,  it  simulates  calendar  time  between 
removals  as  exponential  random  variables.  The  simulation  program  submitted 
with  this  report  may  be  adapted  to  run  as  a FORTRAN  subroutine  of  i.he  JEMS 
program  or  it  may  be  worthwhile  to  rewrite  both  programs  in  a language  more 
suitable  for  representation  of  the  replacement  and  repair  processes  with 
their  different  time  scales,  operating  time  and  calendar  time.  GASP  IV 
is  a candidate  language  (Pritsker  [ 15  ])  since  it  is  a FORTRAN  based 
language  which  can  handle  events  occurring  in  continuous  time  such  as  removals 
in  operating  time  and  events  occurring  in  discrete  time  such  as  repairs  or 
rebuilds  which  need  not  be  monitored  as  accurately  as  is  possible  with  the 
next  event  simulation  method. 
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One  objective  of  the  actuarial  method  is  to  obtain  an  estimate  of 
replacement  requirements.  Obtaining  a good  estimate  of  the  actuarial 


failure  rates  is  only  part  of  the  job.  The  replacement  age  estimate  must 
then  be  transformed  into  an  estimate  of  replacement  requirements.  The  only 

• theoretical  property  true  of  replacement  requirement  estimates  obtained 
from  actuarial  failure  rates,  the  product  limit  estimator,  and  the 
maximum  likelihood  estimator  is  that  they  are  all  maximum  likelihood  esti- 
mators for  their  respective  sample  information  and  replacement  age  models. 
Further  comparison  must  be  done  by  comparing  their  information  matrices  or 
by  empirical  means  such  as  simulating  a replacement  process  and  comparing 
the  results  where  the  engine  life  data  is  generated  from  the  estimators  of 
replacement  age  being  compared.  That  will  be  done  now,  to  illustrate  the 
next  event  method  of  simulation  and  to  compare  the  estimates  of  replacement 
requirements  obtained  from  estimates  of  the  removal  age  failure  rates  and 
cumulative  distribution  function. 

In  order  to  illustrate  the  next  event  method  and  to  compare  it  with 
the  time  slice  method,  20  repeat  runs  of  a 300  hour  flying  program  for  one 
single  engine  aircraft  will  be  simulated  by  hand.  The  number  of  engines 
required  to  achieve  the  required  flying  hours  will  be  collected  for  each  run 
from  which  the  mean,  variance,  percentiles,  and  the  cumulative  distribution 
function  will  be  estimated. 

In  Table  5 engine  ages  at  removal  are  simulated  from  each  of  the 
estimators  including  the  "actuarial  cumulative  distribution  function" 
obtained  from  the  actuarial  failure  rates.  This  function  was  graphed  in 
Figure  3,  Section  5,  part  I.  Table  6 shows  the  calculation  of  engine 
requirements  based  on  simulated  operating  times  from  Table  5.  Figure  15 
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illustrates  the  results  simulated  from  the  theoretical  and  product  limit 
cumulative  distribution  functions.  Table  7 summarizes  the  frequency  data 
from  Table  6 and  graphs  for  comparison  are  in  Figures  16,17,  and  18-  The  product 
limit  estimates  were  as  close  or  closer  to  the  theoretical  values  of  all 
measures  as  any  estimates.  The  actuarial  failure  rate  estimates  didn't 
stand  a chance  because  the  engine  lives  generated  from  actuarial  failure 
rates  used  many  random  numbers  not  used  in  the  other  simulations.  Because  of 
this,  the  actuarial  failure  rates  were  transformed  into  an  actuarial  cumulative 
distribution  function  from  which  engine  lives  were  generated  from  the  same 
random  numbers  used  to  generate  engine  lives  from  the  product  limit  and  theoretical 
distributions.  (The  maximum  likelihood  estimator  mentioned  in  part  I and 
derived  in  Appendix  A is  equivalent  to  the  product  limit  estimator  for 
simulation  of  replacement  requirements.) 

In  Table  5,  Column  1 is  the  index  number  of  the  engines  needed  to  operate 
twenty  planes  over  a period  of  three  hundred  flyinq  hours  each.  Only  ten 
engines  are  shown.  Approximately  65  engines  were  required.  Column  2 
are  the  random  numbers  taken  from  Mihram  [ 2 ].  The  first  random  number  is 
used  to  generate  the  engine  lives  for  columns  3,  4,  and  5.  All  the  random 
numbers  listed  were  used  for  column  6,  the  actuarial  failure  rate  engine 
life.  Column  3 is  the  engines  life  generated  from  the  cumulative  distribu- 
tion function  obtained  from  the  actuarial  failure  rates.  That  cumulative 
distribution  function  is  plotted  in  Figure  3 . Column  4 is  the  engine  life 

generated  from  the  product  limit  cumulative  distribution  function  in 
Figure  4 . Column  5 is  the  engine  life  generated  from  the  theoretical 
cumulative  distribution  function  in  Figure  3 . The  procedure  for  columns 
3,  4,  and  5 is  to  take  the  random  number  in  column  2,  find  it  on  the  vertical 


axis  of  the  graphs  in  Figures  3 and  4,  move  across  until  it  meets  the 
function,  look  down  to  the  time  axis,  and  record  the  removal  time. 

Column  6 is  the  actuarial  failure  rate  removal  time.  These  times  were 
generated  from  Table  2,  the  actuarial  interval  failure  rates.  The  first 
random  number  listed  in  column  2 is  compared  with  the  first  interval 
failure  rate  in  Table  2,  Column  8.  If  the  random  number  is  less  than  the 
interval  failure  rate,  the  actuarial  failure  rate  life  becomes  the  mid- 
point of  the  actuarial  age  interval.  If  the  random  number  is  greater  than 
the  interval  failure  rate,  take  the  next  random  number  in  column  2 and 
compare  it  to  the  second  interval  failure  rate  in  Table  2.  Continue  in 
this  way  until  the  engine  life  time  for  engine  number  one  is  found,  then 
repeat  this  procedure  for  the  next  engine.  For  maximum  time,  the  failure 
rate  is  infinity,  so  if  the  engine  flies  this  long  it  will  have  flown  two 
hundred  hours. 

The  number  of  engines  required  for  twenty  aircraft  varies  slightly  for 
each  method.  Each  of  the  engine  lives  generated  has  been  assumed  to  be  for 
a new  engine,  so  each  time  an  engine  is  removed  a new  engine  replaces  it. 
Table  6,  Column  1 gives  the  number  of  the  aircraft  that  was  used.  Columns 
2,  4,  6,  and  8 are  the  engines  lives  used  by  each  plane  to  fly  three  hundred 
hours  for  the  actuarial  cumulative  distribution  function,  product  limit, 
theoretical  cumulative  distribution  function,  and  actuarial  failure  rates 
respectively.  Columns  3,  5,  7,  arid  9 are  the  numbers  of  engines  required 
for  each  plane  for  the  actuarial  cumulative  distribution  function,  product 
limit,  theoretical  cumulative  distribution  function,  and  actuarial  failure 
rates  respectively.  The  engine  lives  for  each  plane  to  fly  three  hundred 
hours  are  taken  from  Table  5.  The  engine  lives  were  taken  in  chronological 
order  from  Table  5 until  the  total  flying  time  was  greater  than  or  equal 
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to  three  hundred  hours.  For  the  actuarial  cumulative  distribution  function, 
the  first  aircraft  took  six  enqines  to  reach  three  hundred  flying  hours.  The 
first  six  removal  times  from  Table  6,  Column  2 are  33,  90,  47,  33,  16,  and 
200  which  totals  419  hours.  The  excess  time  was  not  carried  over  to  the 
next  aircraft,  so  the  engine  lives  for  the  second  aircraft  are  200  and  200. 

In  Column  3 the  number  of  engines  required  is  the  number  of  engine  lives 
needed  to  reach  or  exceed  three  hundred  flying  hours,  in  this  case,  six 
engines  are  needed  for  the  first  aircraft. 

Figure  15  is  a chart  of  the  product  limit  engine  lives  and  the  theoretical 
engine  lives.  It  shows  when  the  engine  replacements  occurred  over  the  time 
period  of  three  hundred  hours.  Frequently  data  is  tabulated  in  Table  7. 

Column  1 is  the  number  of  engines  that  are  required  for  an  aircraft.  Columns 
2,  3,  4,  and  5 are  the  number  of  aircraft  that  required  the  number  of  engines 
specified  in  Column  1.  Figures  16,  17,  and  18  are  the  cumulative  distribution 
functions  taken  from  the  frequency  table  for  each  of  the  three  ways  engine 
requirements  were  generated.  They  are  compared  to  the  theoretical  engine 
requirements  cumulative  distribution  function. 

A computer  program  has  been  written  for  comparing  estimators  on  the  basis 
of  how  well  they  simulate  engine  requirements  for  one  single  engine  aircraft. 
The  program  does  not  yet  simulate  engine  requirements  from  the  theoretical 
distribution  of  engine  lives  but  this  will  be  incorporated  soon  and  full 
scale  testing  could  be  done  for  larger  samples  than  used  in  this  section. 


TABLE  5 


RAN DO! 1 ENGINE  LIVES  FOR  REPLACEMENT  SIMULATION 


Col umn 
Number 


Col umn 
Name 


1 

2 

3 

4 

5 

C. 

Engine 

Number 

Random 

Number 

Actuarial  c.d.f. 
Engine  Life 

Product  Limit 
Li  fe 

Theoretical 
Li  fe 

Actuarial 
Failure  Rate 
Life 

1 

.153 

3.3 

10 

38 

20 

2 

.337,-079 

.165 

90 

92 

85 

100 

3 

.203 

47 

33 

52 

20 

4 

.162 

33 

33 

41 

20 

5 

.078 

16 

10 

21 

20 

6 

.876, .58, 
.09 

200 

200 

200 

100 

7 

. //8, . 30, .46, 
.35,  .10 

200 

200 

200 

180 

8 

.909, .79, .45, 
.69, .18 

200 

200 

200 

180 

9 

.931, .63, .53, 
.10,-43 

200 

200 

200 

200 

10 

.326, .92, .56, 

89 

82 

82 

180 

.63, .07 


53 


200+100 


TABLE  7 


FREQUENCY  DATA  ON  ENGINE  REQUIREMENTS 


Col  umn 


Number 

1 

2 

3 

4 

5 

Col umn 
Name 

Number  of 
Engines 
Required 

Actuarial 

c.d.f. 

Frequency 

Product 

Limit 

Frequency 

Theoretical 

c.d.f. 

Frequency 

Actuari al 
Failure  Rates 
Frequency 

2 

8 (.4) 

7 

(.35) 

7 (.35) 

9 (.45) 

3 

6 (.3) 

7 

(.35) 

7 (.35) 

9 (.45) 

4 

4 (.2) 

1 

(.05) 

3 (.15) 

1 (.05) 

5 

0 

2 

(.1) 

0 

0 

6 

2 (.1) 

3 

(.15) 

2 (.1) 

0 

7 

_0 

_0 

_L  (.05) 

_1  (-05) 

20 

20 

20 

20 

90th  Percentile 

6 

6 

7 

7 

80th  Percentile 

6 

6 

6 

4 

70th  Percentile 
Mean 

4 

3.10 

6 

3.35 

6 

3.30 

3 

2.75 

Variance 

1.57 

2.13 

2.22 

.93 

Standard 

Deviation 

1 .25 

1 .46 

1 .49 

.97 
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9.  The  Actuarial  Poisson  Approximation  to  Spare  Engine 
Requirements  and  Other  Analytical  Approximations 

It  is  not  sufficient  for  engine  management  to  predict  only  the  expected 
number  of  engine  replacements  that  require  repair  or  rebuilding.  The  Air 
Force  must  plan  spares  requirements  to  give  a comfortably  high  probability 
of  meeting  mission  requirements,  usually  80'  or  90",  i.e.,  the  number  of 
spares  must  be  sufficient  to  meet  replacement  requirements  for  the  flying 
hour  program  with  probability  .8  or  .9.  In  order  to  determine  these 
numbers,  the  80th  and  90th  precentiles  of  the  replacement  requi rements 1 
distribution  must  be  estimated.  The  actuarial  method  contains  a procedure 
to  do  this.  The  simulation  program  can  be  easily  revised  to  produce  these 
estimates.  Other-  analytical  approximations  may  be  possible,  perhaps  com- 
bined with  simulation.  Each  of  these  approaches  will  be  described  in  this 
section  and  the  foundations  for  the  Poisson  approximation  will  be  reviewed 
to  show  a potential  improvement  in  its  accuracy. 

The  actuarial  method  for  estimating  safety  stock  levels  in  AFM  400-1 
Chapter  8 [ 13  ] is  based  on  the  assumption  that  demand  for  spare  engines 
has  a Poisson  distribution  with  rate  parameter  or  expected  demand  per  unit 
calendar  time  equal  to  the  rate  that  engines  require  repair  or  rebuilding. 

(If  there  are  fewer  spares  than  will  eventually  be  in  repair  or  are  being 
rebuilt,  then  some  aircraft  will  not  have  engines.)  The  Safety  Level  Table 
Figure  8-1  of  AFM  400-1  is  a table  of  the  90th  percentile  of  a Poisson  dis- 
tribution for  different  values  of  its  parameter.  The  Poisson  assumption  is 
unlikely  to  be  true  exactly,  but  the  conditions  for  use  of  a Poisson  approxi- 
mation may  be  near  realization.  Those  conditions  are  (Cox  [ 12  ] Chapter  6): 
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1.  There  are  many  simultaneously  operating  or  installed  engines,  and 

2.  The  distribution  of  the  residual  calendar  times  from  engine 
installation  to  removal  is  the  same  for  all  engines  and  those  times 
are  independent. 

Condition  2 may  be  violated  for  multiple  engine  planes  and  it  may  also  be 
violated  because  some  engines  are  not  new  at  installation. 

The  conditions  above  should  be  tested  if  the  Poisson  approximation 
for  spare  engine  requirements  is  to  be  used.  If  the  conditions  are 
realized,  then  the  Po  jn  approximation  should  be  modified  to  take  into 
account  the  initial  failure  rate  after  installation.  Unless  the  calendar 
time  between  engine  installation  and  removal  is  exponentially  distributed, 
the  initial  failure  rate  can  be  used  to  improve  the  Poisson  approximation. 
(If  the  initial  failure  rate  is  higher  than  average,  more  spares  are 
required.  ) 

Meanwhile,  the  next  event  simulation  of  replacement  requirements  can 
give  all  information  about  the  distribution  of  replacements  including 
percentiles.  For  instance  the  90th  percentile  of  simulated  replacement 
requirements  are  obtainable  from  Figures  16  , 17,  and  18  , and  they  are 

shown  in  Table  7.  They  are  obtained  by  reading  the  figures  across  from 
the  vertical  axis  at  probability  value  0.90  and  reading  down  to  the 
horizontal  axis  to  obtain  the  value  of  engine  requirements  that  was  not 
exceeded  in  90  of  the  simulation  runs.  The  accuracy  of  this  estimate  could 
be  estimated  from  many  repeated  simulation  runs. 


Another  approximation  to  the  number  of  replacement  requirement,  is 
possible  by  an  adaptation  of  the  central  limit  theorem.  It  is  known  that 
the  number  of  replacements  with  new,  independent,  identical  components  has 
approximately  a normal  distribution  after  a long  time,  regardless  of  the 
life  distribution  of  the  components  (Cox  Chapter  3 [ 12  ]).  The  central 
limit  theorem  is  remarkably  robust  and  will  tolerate  some  variation  in  the 
residual  operating  time  distribution  of  the  components,  (Lindberg  Conditions, 
Gnedenko,  Chapter  8 [ 23  ]).  This  tolerance  is  probably  sufficient  so 
that  the  normal  approximation  to  replacement  requirements  is  still  ade- 
quate for  engines  even  though  some  engines  are  not  new  at  installation. 

This  long  time  approximation  could  be  combined  with  a simulation  to 
determine  replacement  requirements  for  the  near  future.  Some  additional 
research  is  required  here  to  develop  the  normal  approximation. 

The  temptation  is  always  present  to  assume  a simple  model  of  engine 
operation  for  convenience  in  predicting  replacement  requirements.  It 
would  also  be  desirable  to  incorporate  additional  information  about  an 
engine's  history  into  a model  of  engine  operation  to  obtain  more  accurate 
prediction  of  removal  times  and  replacement  requi rements.  Several 
approaches  are  possible  that  should  be  explored. 

1)  Postulate  that  the  failure  rate  function  of  the  rer  r il 
times  (regardless  of  whether  the  removal  is  for  repair 
or  rebuild)  is  of  the  form  (Cox  [ 14  ]) 
r(t)  = rQ(t)e~? 

where  Z is  a row  vector  of  concomittant  variables  with  some 
relation  to  engine  age  at  removal  and  is  a column  vector 
of  constants.  Likely  concomittant  variables  are  the  number 
of  prior  repairs  and  the  age  at  last  installation. 
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2)  Postulate  that  the  process  causing  combined  removals  interacts 
with  an  age  dependent  probability  that  determines  whether 

the  engine  is  repaired  or  rebuilt.  (Older  engines  are  more 
likely  to  be  rebuilt  than  nearly  new  engines  when  removed.) 

3)  Postulate  that  the  replacement  requirements  process  is  a 
Markov  process  in  operating  time  and  that  there  is  a known 
random  transformation  from  operating  time  to  calendar  time 
which  may  depend  on  the  Markov  process.  The  number  of  re- 
placements required  is  also  a Markov  renewal  process  under 
rather  general  conditions. 

Each  of  these  three  postulated  models  has  properties  that  make  them 
useful  for  different  aspects  of  engine  management.  The  first  model  is  con- 
venient to  estimate  and  can  be  used  to  determine  trends  in  the  performance 
of  engines  (Tarone  [18]).  The  second  model  is  convenient  for  determining 
operating  time  between  engine  overhauls.  The  third  model  converts  replace- 
ments occurring  in  operating  time  to  replacements  in  calendar  time  as  a 
function  of  the  flying  hour  program  and  the  rate  at  which  each  aircraft 
contributes  to  the  program. 


SECTION  3 


RECOMMENDATIONS 

Recommendation  1:  Do  not  smooth  the  "crude"  actuarial  failure  rates. 

The  smoothing  procedure  is  a vestige  of  the  life  insurance  industry 
where  there  was  no  reason  to  assume  human  failure  rates  were  not  smooth. 

The  converse  is  true  for  engine  failure  rates  due  to  component  failure  and 
the  high  probability  of  removal  at  inspection.  To  smooth  the  failure 
rates  is  to  suppress  the  information  most  useful  for  engine  diagnosis  and 
for  prediction  of  replacement  requirements.  Recommendation  1 can  be  adopted 
immediately  and  will  only  reduce  computation. 

The  actuarial  method  employs  a linear  extrapolation  of  earlier  failure 
rates  where  data  grows  sparse  and  the  actual  data  is  not  used.  This  data  can 
and  should  be  used  in  calculation  of  the  engine  removal  time ''distribution 
and  failure  rates.  Several  recommendations  to  follow  will  use  all  data. 

Even  if  none  of  them  is  adopted,  it  is  still  possible  to  get  enough  observa- 
tions in  each  actuarial  age  interval  to  estimate  meaningful  failure  rates  by 
making  some  of  the  intervals  larger.  The  ultimate  extension  of  variable 
actuarial  age  interval  estimators  is  the  estimator  with  exactly  one  removal 
per  interval.  Because  this  estimator, call ed  the  "product  limit"  estimator, 
was  already  available,  no  further  study  of  actuarial  type  estimators 
was  done.  Unfortunately  this  estimator  does  not  take  into  account  the  age  of 
engines  at  the  beginning  of  some  calendar  period  of  operation,  information 
which  now  may  be  useful  for  some  purposes.  This  leads  to  two  recommendations 
for  further  development. 

Recommendation  2:  Develop  actuarial  estimators  with  flexible  age 

intervals  and  modify  the  product  limit  estimator  to  take  into  account 
engine  ages  at  the  beginning  of  some  period  or  their  ages  at  the  latest 
installation  if  they  were  installed  after  the  beginning  of  the  period. 
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Recommendation  3:  Determine  whether  it  is  worthwhile  to  use  the 

additional  information  on  ages  for  all  actuarial  purposes. 

The  estimator  with  one  removal  per  age  interval  uses  all  engine  operating 
time  information,  whether  a removal  time  or  the  age  of  a still  operating 
engine, to  obtain  the  cumulative  failure  rates  or  the  cumulative  distribution 
function  of  engine  removal  times.  This  type  of  estimator  also  can  provide 
every  information  product  now  computed  from  the  actuarial  failure  rates, 
and  the  resulting  products  will  usually  be  more  accurate  because  the  actual 
removal  time  information  is  fully  used.  The  product  limit  estimator  requires 
more  computation  than  the  actuarial  estimator,  but  the  required  computation 
is  well  within  the  capability  of  current  computers. 

The  second  recommendation  requires  applied  statistical  research  and  the 
third  requires  some  testing.  If  the  second  recommendation  is  successful, 
there  is  no  reason  to  pursue  the  third.  There  is  every  reason  to  expect 
success  on  the  second  recommendation  because  there  are  possible  estimators 
that  exist  or  could  easily  be  developed  to  use  all  the  information  now 
available  on  engine  ages  and  removal  times.  As  part  of  this  research  contract, 
the  maximum  likelihood  estimator  for  an  engine  removal  time  model  that 
distinguishes  usage,  inspection  and  max  time  removals  was  derived  for  a 
sample  of  removal  times  and  survivors  ages.  The  derivation  can  be  modified 
to  incorporate  engine  ages  at  the  beginning  of  some  calendar  period  or  at 
last  installation.  That  is  recommendation  4. 

Recommendation  4:  Derive  the  maximum  likelihood  estimator  for  the 

removal  time  model  with  statistically  independent  usage  and  inspection 
removal  times.  Test  whether  this  model  is  appropriate  for  engine 
removals . 
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If  the  model  is  appropriate,  an  estimator  will  be  available  that  is 
expected  to  be  simpler  to  use  than  the  product  limit  estimator,  and  it  also 
gives  an  estimate  of  usage  removal  times,  useful  for  engine  diagnostics 
and  for  curve  fitting  attempts  to  identify  a parametric  model  of  engine  usage 
removal  times. 

One  question  that  should  be  addressed  by  the  actuarial  statistical  test 
is  - does  the  official  failure  rate  reflect  the  performance  of  the  current 
population  of  engines?  The  answer  is  provided  by  comparing  the  data  set 
from  which  the  official  failure  rates  were  computed  with  the  data  set  con- 
taining current  engine  removal  data.  The  updating  program  submitted  with 
this  report  provides  a statistically  valid  measure  that  the  two  data  sets 
really  came  from  the  same  population.  It  is  designed  to  compare  all  past 
data  sets  with  current  data  on  removal  times  and  survivors'  ages. 

Recommen..  "ion  5:  Adopt  the  engine  removal  data  updating  program 

submitted  with  this  report  to  replace  the  current  statistical  test  for 
changing  the  official  failure  rates  age  interval  by  age  interval. 

This  updating  program  does  not  use  ages  of  engines  at  the  beginning  of 
any  period  or  at  installation.  A modification  has  been  proposed  that  will 
do  this  but  it  requires  further  development. 

Recommendation  6:  Develop  the  modified  statistical  two  sample  Wilcoxon 

rank  test  to  incorporate  engine  removal  times,  survivors'  ages  at  the 
end  of  the  operating  periods,  and  ages  at  the  beginning  of  the  period  or 
at  installation  whichever  is  later. 

The  actuarial  procedure  for  computing  confidence  limits  for  spare  engine 
requirements  recognizes  that  demands  do  not  always  equal  expected  demand,  so 
additional  spares  are  required  to  keep  the  probability  of  running  out  to  an 


acceptable  level.  However  the  procedure  is  based  on  the  assumption  that 
demands  have  a Poisson  distribution.  This  is  true  only  under  very  special 
conditions  or  when  the  number  of  operating  engines  is  very  large  and  all 
replacements  are  new.  The  Poisson  assumption  should  be  revaluated  with 
two  possibilities  in  mind.  First,  the  Poisson  approximation  may  be  strengthened 
by  a second  order  term  involving  the  initial  failure  rate.  Second,  the 
normal  approximation  may  be  adapted  to  provide  a better  estimate  of  engine 
requirements  after  the  engines  have  been  in  use  for  some  time.  These 
observations  lead  to  recommendation  7. 

Recommendation  7:  Revaluate  the  spare  engine  requirements  computation 

now  based  on  the  assumption  of  a Poisson  demand  distribution. 

Other  data  related  to  engine  removal  times  is  available  and  should  be 
used  to  obtain  more  accurate  estimates  of  engine  removal  times  and 
consequently  of  replacement  requirements;  data  such  as  sorties,  engine  cycles, 
prior  repairs.  X-ray  and  chemical  analyses.  Development  and  testing  of  more 
compreriensive  statistical  estimators  is  necessary  for  significant  improvement 
beyond  the  expected  improvement  from  recommendations  1-6. 

Recommendation  8:  Support  the  research  necessary  to  develop  more 

comprehensive  estimators  of  engine  removal  times. 

Presently,  simulation  is  the  only  statistically  valid  way  to  estimate 
engine  requirements.  The  simulation  program  now  used  to  predict  engine 
requirements  simulates  each  engine  through  each  actuarial  age  interval. 

The  result  is  a time  consuming  program  which  gives  only  an  estimate  of 
expected  replacement  requirements  but  no  confidence  limits  and  no  indication 
of  the  accuracy  of  the  requirements  estimate.  The  next  event  method  can 
simulate  fleet  operation  many  times  and  obtain  information  on  confidence 
limits,  accuracy  of  estimates,  and  even  estimates  of  when  engine  removals 
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will  occur.  The  method  can  be  used  whether  actuarial  failure  rates 
used  or  the  product  limit  estimator  of  recommendation  2 is  used. 

Reconmendation  9:  Change  the  simulation  program  for  replacement 

requirements  from  an  interval  by  interval  simulation  to  a next  event 
type  simulation. 

The  program  prepared  for  this  report  contains  the  ingredients  fu. 
accurate  simulation  of  engine  removals  but  does  not  make  any  pretense  of 
accurately  representing  the  repair  and  rebuilding  of  engines. 

This  has  already  been  done  in  the  JEMS  program  developed  at  AFLC. 

Recommendation  10:  Combine  the  next  event  type  simulation  of  the 

engine  removals  and  the  JEMS  program  simulating  repair  and  rebuilding. 
This  will  have  to  be  done  in  close  cooperation  with  AFLC  because  of  their 
thorough  knowledge  of  the  repair  process  for  engines. 

The  resulting  program  is  likely  to  be  the  largest  program  that  can  be 
handled  in  a reasonable  amount  of  computer  time.  There  is  still  a need  for 
more  accurate  and  more  convenient  methods  to  predict  replacement  requirements. 
A normal  approximation  conditioned  on  the  ages  of  the  present  stock  of 
engines  should  be  developed  for  estimating  engine  requirements. 

Recommendation  11:  Support  the  research  and  development  necessary 

for  a simulation  of  short  term  engine  requirements  combined  with  a 
long  term  estimate  based  on  a normal  approximation. 
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SECTION  I 


r 


INTRODUCTION 


The  ultimate  objective  is  to  derive  the  maximum  likelihood  estimator 
of  the  cumulative  distribution  function  (c.d.f.)  of  engine  removal  times 
from  a progressively  censored  sample.  The  c.d.f.  is  assumed  to  be  of 
the  form 

) \ (1) 


_ ( i(t) 

F(t)  = i-F1(t)  J n (i-Pi 

( i-1 


where  F^(t)  is  the  tail  c.d.f.  of  another  non-singular  c.d.f.  possibly 
truncated  and  representing  usage  removals,  and  the  p^  are  the  unconditional 
probabilities  of  removal  at  the  i*^  inspection  time  (fixed)  where  i(t)  is 
the  index  of  last  inspection  prior  to  or  at  time  t.  Throughout  this  deriva- 
tion time  is  measured  in  flying  hours,  not  calendar  time.  Progressive 
censoring  occurs  because  estimation  of  the  c.d.f.  takes  place  at  fixed 
calendar  times,  and  in  addition  to  data  on  the  flying  hours  at  removal 
of  engines  that  have  been  replaced,  the  current  accumulated  flying  hours 
of  installed  engines  that  have  not  yet  been  removed  is  also  known. 

The  plan  of  this  derivation  is  in  three  stages.  First,  for  the 
sake  of  review,  the  empirical  distribution  function  for  an  ordered  random 
sample  of  engine  removal  times  T^  < T2  T^ 


i 

i 

* 
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will  be  derived  by  the  method  of  maximum  likelihood.  Then  the  maximum 
likelihood  estimators  of  F^(t)  and  the  Rj , j = 1,  2,...,  k,  will  be 
derived  for  an  uncensored  sample  of  engine  removal  times.  Substitution 
of  those  estimators  into  the  c.d.f.  (1)  yields  the  empirical  distribution 
function  as  an  e.otimate  of  the  model  (1)  but  with  additional  information 
about  usage  removals  and  inspection  removals.  Last,  the  maximum  likelihood 
estimators  of  F^(t)  and  the  pj are  derived  from  a progressively  censored 
sample.  Substitution  into  (1)  yields  the  "product  limit"  estimator  of 
Kaplan  and  Meier  [l]. 


SECTION  LI 


MAXIMUM  LIKELIHOOD  DERIVATION  OF 


THE  EMPIRICAL  DISTRIBUTION  FUNCTION 


Given  an  ordered  random  sample  from  an  arbitrary  c.d.f.,  it  is 
frequently  assumed  that  the  mass  of  the  estimator  will  be  placed  only 
at  observed  data  in  the  sample.  Then  the  likelihood  function  of  the 
sample  is  the  product  of  masses  f , i=l,2,...,N,  at  the  data  values 

t1  — C2  ~ ' ' ' — tN’ 


L(ti,...,t  i F(t))  = n f 
w i=l 

and  the  masses  are  subject  to  the  constraints  of  non-negativity  ana 


I f . = 1 


so  that  the  estimator  of  F(t)  is  a c.d.f.  The  values  of  f that 
maximize  the  likelihood  function  are  all  equal  to  l/N.  This  can  be 


verified  by  the  Lagrange  multiplier  method  or  by  substitution  of  the 


constraint 


f.T  = 1 - £ f. 

N i=l  1 

into  the  likelihood  function  and  setting  derivatives  of  the  log  likelihood 
function  with  respect  to  f ^ equal  to  zero. 


6 log  L _ 6_  { £ log  f . + log  (1-  £ f ,)  } 

6f ± 6f  ± j =1  J j =1  J 


k 


SECTION  III 

MAXIMUM  LIKELIHOOD  DERIVATION 
OF  THE  FULL  SAMPLE  ESTIMATOR 

The  engine  removal  times  Tl’ • ' ’ >TN  in  the  sample  are  assumed  to 
be  independent  and  identically  distributed  according  to  F(t) , (1).  In 
addition,  the  engine  removal  code  specifies  whether  the  removal  was 
during  usage,  at  inspection,  or  at  maximum  time.  Define  the  following 
sequences : 


(nj),  the  number  of  removals  at  the  jth  inspection  and 
tmax>  j = l>2,. . . ,k+l;  (t^^  = tfc+1) 

{m.},  the  number  of  usage  removals  between  the  (j-l)j^t  and  the 
jth  inspection; 

{M  },  the  cumulative  number  of  usage  removals, 

i 

M.  = I m. , j=l,2, . . . ,k+l ; 

J i=l  1 

{ t ^ } , the  sequence  of  usage  removal  times,  i=l , 2 , . . . ,Mk+2 ; and 

f 

{ t . } , inspection  times  and  t 
J max 

These  sequences  turn  out  to  be  sufficient  statistics.  The  sample 
size  N can  be  expressed  in  terms  of  usage  removals,  inspection 
removals  and  maximum  times  removals  as 


k+1 

N = Mk+1  + I n 

j = l 


The  likelihood  function  for  the  sample  information  is 


\+ir  J(ti}  1 k r J_1  1 ni 

n f,(t.)  n (l-p.)  n p n (i-P.)  (i-f  , (t ' ) ) J 

i=i  L 1 j=i  J J j=i  L J i=i  1 1 j J 


[h-hUnax  > ^ Vl 


with  f^(t^)  denoting  the  density  of  F^(t),  discrete  or  continuous,  at  t^, 

t “ = t -e  for  small  e>0,  and  the  notation  i(t.)  is  used  to 
max  max  l 

indicate  the  index  of  the  last  inspection  prior  to  the  usage  removal 
at  t^.  The  first  term  of  the  likelihood  function  is  the  probability 
of  the  usage  removals,  the  second  term  is  the  probability  of  inspection 
removals,  and  the  last  is  the  probability  of  all  the  survivals  to 


t . C is  a combinatorial  constant, 
max 


In  order  to  make  the  likelihood  function  positive,  all  of  the 
f^(t^)  should  be  positive,  and  to  maximize  it,  l-F^(tj)  and 
1-Fj_(t  ) should  be  as  large  as  possible  consistent  with  the  constraints 

on  the  distribution  function  F(t);  F(0  )=0,  F(“>)=1,  and  F(t)  non- 
decreasing in  t.  As  in  the  maximum  likelihood  derivation  of  the 
empirical  distribution  function,  all  the  mass  of  F^(-)  should  be  placed 
at  the  observations  {t,}  and  t . 

1 ItlaX 


Then 


M 

l-FL(t;)  = 1 - f1(t.) 

i=l 
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and 


1-F  ( t ) = 1 - 
l max 


However,  not  all  of  the  jumps  in  F^(')  will  be  of  the  same  size  as  in 
the  empirical  distribution  function.  The  derivative  of  the  natural 
logarithm  of  the  likelihood  function  with  respect  to  f-^(t^)  gives 


M. 

6 log  L/ifrti)  = l/f1(ti)  - f1(ti)/  U - f^t^  > 

1 1 i=l 

for  all  i=Mj_i  + 1,..,M^,  j=l,2, . . . ,k+l.  This  indicates  that  the 

maximum  likelihood  estimators  of  the  jumps  will  be  the  same  in 

intervals  before  the  first  inspection,  between  inspections,  and  after 

the  last  inspection.  These  jump  sizes  will  be  denoted  by  ■ * • *ifc+l* 

The  likelihood  function  now  simplifies  to 

k+I  r i-1  -|  11  ~ H k+1  r j-1  . j n 

n \ ii  n (i-Pi)  i i-i  non  (i-Pi){i-  r ^(m.-m^i) 

i=iL  j-i  3 J j=i  L J i=i  i-i  J 

where  = 1-  The  log  likelihood  function  is 

k+1  i-1 

Z (Mi-Mi_1){ln  ft  + l ln(l-p . ) } + 
i=l  j=l  3 

k+1  j-1  j 

+ Z n.  { In  p.  + Z ln(l-p.)  + ln(l-  £(M.-M  .)f^)  ) + In  C 

J-1  J 3 i-1  i“l  1 

Take  the  derivative  of  the  log  likelihood  function  with  respect 
to  p.  and  set  it  equal  to  zero,  j=l,2,...,k. 
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k+1 


•1 


n . 

J 


\SSJLT  E (M.-Mi_i)  (j-~)  + - 

6 P-i  i=j+l  J J 


k+1  i 

2 t 


= 0 


Remove  the  1-p . terms  from  the  summations  to  obtain 
J 


"j  ! k+1 

p7  = T^T  ,2.+1  (Mi~Mi-l  + ni} 
J 1=3+1 


or 


k+1 

nj  = pj  { Mk+i  - Mj  + ^ ni  } 

This  gives  the  maximum  likelihood  estimator 

k+1 

Pj  = V {Mk+l  " Mj  + ,£.  ni} 

i=J 

which  may  be  rewritten  as 

3“1 

Pi  = n.J  {N-M,-  E n.}  j=l,2,...,k 

3 J i=l 

where  the  denominator  is  the  number  of  survivors  to  the  jth  inspection. 

Thus  pj  is  simply  the  proportion  of  survivors  removed  at  the  jUi  inspection. 

Take  the  derivative  of  the  log  likelihood  function  with  respect 
to  fj  and  set  it  equal  to  zero,  i=l , 2 , . . . ,k+l , 


6 In  L 1 k+1  n (M.-M  ) 

=(M  -M  ) — - Y h 1 1-1 

Sf.  1 i i-Tf.  uL 


h=i 


h 

Z 

3=1 


1-  E 


Divide  out  to  obtain  the  general  relation  for  the  maximum 

likelihood  estimators  of  f^ 
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I 


k+1 

— = l 

fi  h=i 


II 

1-  l (M.-M.  ,)f. 
j = l J ^ J 


1-1,2,. 


Suppose  for  the  moment  that  there  is  only  one  inspection  time,  k=l. 

Then  to  obtain  the  maximum  likelihood  estimators  of  f^  and  f , the  equations 


fi  ^Vi 


1 m1f1-m2f2 


f2  1-,nifrra2f2 


must  be  solved  simultaneously.  The  first  equation  contains  l/f2  as  the  last 


term,  and  rearrangement  gives 


1~mifrnifi  i 

fl(1-lfl)  = f2 


The  second  equation  can  be  solved  for 


f 2~fA 

2 n2+m2 


Combining  these  two  expressions  for  f2  gives 
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fj  = 1/ (n1-»-n2-Hn1-Kn2 ) = 1/N 


f2  = (l/N)  (n1-*-n2-»-in2)  / (n2+m2)  = (1/N)  (l-p1> 


Now  return  to  the  general  problem  of  solving  the  system  of  k+1  equations 
(3)  simultaneously.  The  answer  is 


f,  = (1/N)  { n (l-p.)}' 
J-l  3 


The  equation  (3)  yields  a recursion  relation 


% (YVl>£J 


+ 4-1 


which  may  be  used  to  verify  the  estimators  (4) . This  has  been  done  for 
k=2  and  is  very  tedious,  and  the  result  (4)  for  larger  k will  be  verified  by 
induction  on  i for  an  arbitrary  number  of  inspections  k>i  from  recursion 
relation  (5) . 

If  it  can  be  shown  that  recursion  relation  (5)  yields 


7-7  = Np  n (1-p  ), 

i i+1  J-l  11 


the  induction  is  complete  because  the  difference  above  is  equivalent  to  the 


1 


solution  (4) . The  rest  of  this  proof  is  to  verify  that  the  center  ratio  of 
(5)  is  equal  to  the  right  hand  side  of  (6) , 


^ i-1 

Np.  n 
1 j-i 


<1-pj> 


1-  \ m.f . 

j-1  3 3 


(7) 


Substitute  the  solution  (4),  j=l,2,...,i,  for  t into  (7).  This  is  fair 
because  under  the  induction  method,  the  result  (4)  is  assumed  true  for 
j=l,2,...,i,  and  is  to  be  shown  for  i+1.  This  substitution  yields 


„ i-1  . n N 

Np  n (1-p  ) = 

J-1  3 i J-1  - 

N-  l [m  / n (1-p  )] 
j-1  3 k=l  k 


1 1-1 
N-  l m,-  l n 

J-l  3 3 

i-1  „ i-1  - 

(N-m  ) n (1-p  )-m  n (1-p 
1 k=l  K k=2 


to  be  verified.  This  will  be  done  term  by  term.  Each  product  in  the  denom- 

i 

inator  includes  1 which  will  have  coefficients  N-  £ m exactly  as  in  the 

j=l  J 

numerator.  It  remains  to  show  that  the  rest  of  the  terms  in  the  denominator 


give  £ n-'  Use  t^e  exPansi°n 


i-1.  i-1. 


11  (1-pk^=1~  ^ pi+  ^ piPi_,,'(_1^  11  pk 

k=l  k k=l  1 1*k  1 3 k-1  K 


to  obtain  the  coefficients  of  p^  in  terms  of  N,  n^,  and  p^  for  jj<k. 


Starting  with  k=l,  the  p..  term  is 


(N-m1>  (-p^-nj. 


The  p2  term  is 


(N-mi)(-p1)+(N-mi)p1p2-m2(-p2)=(N-m1-m2-n1)(-p2)=  -nr 


Similarly  the  p^  term  is 


(N-n^)  [-p3+PlP3+  p2p3-p1p2p3]-m2[-p3+p2p3l  -m3[-p3]=  -n^, 


S 

/ 


and  it  is  hoped  that  this  is  sufficient  demonstration  that  the  coefficients 

A 

of  p^,  k=l,2, . . . ,i-l  when  collected  as  above  yield  with  the 

negative  sign  just  as  in  the  numerator. 

It  is  of  some  interest  to  estimate  the  jump  at  t resulting  from  the 

max 

maximum  likelihood  estimator  f^,  i=l,2, . . . ,k+l.  That  jump  is  estimated  by 


^ - k+1  k+1  m 

Fl(t  max)  = 1_,Lmifi=  1_  N I-I  7 

t_1  1 1 JI  (1-p.) 

j=l  J 


For  k=l,  this  is 


2 . 1 


a very  reasonable  estimator  since  n„/N  is  the  fraction  of  survivors  to  t 

1 2 max 

and  (1-p  )-1  represents  the  conditioning  on  surviving  the  first  and  only 
inspection . 

A similar  result  may  be  proved  for  an  arbitrary  number  inspections  k by 
use  of  the  method  that  verified  the  form  of  f^  (4).  Equation  (9)  becomes 


k k „ k 

' N'1  l'1 1 «-”1)i;i<1-P1>-"2i;2<1-^>-  ' -W 


when  the  product  is  factored  from  every  term  in  the  sum.  Relation  (8)  is 


then  used  to  show  that 
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and  the  conditional  probability  interpretation  again  yields  the  interpretation 
that  this  estimator  is  the  sample  proportion  of  removals  n^^M  conditioned 
on  surviving  all  inspections. 

It  is  reassuring  to  find  that  when  the  estimators  for  F^(t)  and 
p , j*l,2, , . . ,k  are  substituted  into  the  assumed  cdf  (1),  the  resulting 
estimator  of  engine  removal  time  cdf  is  the  empirical  cdf  (2) . This  is 
easily  verified  for  a small  number  of  inspections  k.  For  k=0,  the  result 
is  immediate  because  the  cdf  jumps  by  1/N  at  each  removal  time.  For  k=l, 

i(t) 

F(t)  = 1-F  (t)  n (1-p  ) = 
i=l 

» 

/ i/N  for  t^t<t^+^  an<*  t<t^,i£m^ 


The  third  term  above  simplifies  to  (n^+m^+i)/N,  the  proportion  of  removals 
at  or  prior  to  time  t,  so  F(t)  is  the  empirical  cdf. 


This  result  may  be  proved  for  an  arbitrary  number  of  inspections  k by 
another  application  of  the  relation  (8).  First  verify  that  F(t)  is  h/N  prior 
to  any  inspection  where  h is  the  number  of  usage  removals  that  have  occurred, 


F(t)  = 1-U-hfJ  = l-(l-h/N)  = h/N 
by  simple  substitution  for  f^. 

r 

Then  verify  the  result  for  an  arbitrary  inspection  time  t^; 


= 1-[1-  J m.f.]_n  (1-p.) 


i=l  i=l 


j i“l  J 

= 1-  [1-  l [m  / n (l-pk) ] ] n (1-p.). 

i=l  k=l  i=l 


Factor  the  product  IT  (l~Pk)  from  the  denominator  of  the  F^t)  term  and  cancel 


common  terms  of  the  two  products,  resulting  in 


F(t)  = 1- [ (N-m  ) n (l-p,)-m9  n (1-p  )-...-m  ] . 

1 k=l  k=2  K J N 


Use  relation  (8)  and  the  value  of  p^  to  obtain 


j i 

F(t)  ® 1-  -|[N-  l (m1+ni)]  = l (m^^^), 


I 


(The  denominator  of  1-p^  cancelled  the  numerator  of  F^(t).) 

Last,  verify  that  F(t)  is  the  empirical  cdf  for  a t between  inspections, 

i i 

t .<t<t , 

3 1+1 


F(t)  = 1- [1-  1 m f -hf  ] H (1-p,), 
i=l  1 1 J i=l 


for  h<m.  , and  t <t<t  ,,  , . 

— l+l  m.+h  — m.+h+l 

1 3 


Substitute  for  f.,  factor  out  II  (1-p.)  from  the  denominator  of  the  F (t) 
1 k=l  k 

term,  and  cancel  it  with  the  product  in  the  numerator  to  obtain 


F(t)  = 1-  — [(N-m  ) n (1-p  )-m  H (1-p  )-. 

N 1 k=l  k 2k=2  k 


. ,-m.(l-p.)-h] 


= I I ("'i411! 

i=l 


)+h]/N 


the  empirical  cdf  from  the  relation  (8)  slightly  modified  to  account  for 
the  fact  that  the  product  index  ranges  up  to  j . 

The  appearance  of  the  maximum  likelihood  estimator  of  (1)  from  a full 
sample  is  similar  to  Figure  1 of  the  empirical  c.d.f.  except  that  jump  sizes 
increase  where  there  are  multiple  observations  at  the  same  time  such  as 
occur  at  inspection  times.  It  is  more  interesting  to  note  that  the  esti- 
mator of  F^(t)  is  similar  to  the  empirical  c.d.f.,  but  the  jumps 


95 


vary  in  height  from  one  interinspection  interval  to  the  next  because  fewer 
engines  survive  the  inspection  and  go  on  to  usage  removal  in  the  next  inter- 
inspection interval.  For  example,  suppose  the  usage  removal  times  distributed 
according  to  F^(t)  were  observed  to  be 

{t1>  = {10,  15,  25,  40,  80) 

and  suppose  an  inspection  was  held  at  50  hours  when  two  engines  were  removed 
and  three  more  were  removed  at  max  time, 100  hours.  The  estimator  of  inspec- 
tion removal  probability  is 

P = 2/6  = .333 

and  for  the  estimate  of  F^(t) 

fx  = 1/10  = .10 

and 

f2  = (1/10) (1-. 333)-1  = .150 

The  estimate  of  F^t)  is  shown  in  Figure  2.  The  mass  at  max  time  is  included 
in  the  estimate  of  F^t)  but  that  is  optional. 
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SECTION  IV 


MAXIMUM  LIKELIHOOD  DERIVATION  OF  THE  PROGRESSIVELY 
CENSORED  SAMPLE  ESTIMATOR 

A progressively  censored  sample  contains  values  of  lninlX^,,  Y.  ) , 1=1,2, 

..,N  where  the  X^  are  the  random  variables  whose  c.d.f.  we  wish  to  estimate 
and  the  are  nuisance  variables.  This  situation  can  arise  either  when  all 
N items  begin  testing  at  once  and  some  are  removed  from  testing  before  failure 
at  times  Y^  or,  as  in  Air  Force  engine  management,  the  sample  data  taken  at 
fixed  times  include  some  replacement  times  X^  and  some  ages  Y^  of  engines 
which  have  not  yet  been  replaced.  The  age  Y^  of  the  survivors  provides  ad- 
ditional information  about  the  c.d.f.  of  the  replacement  times  X^  and  should  be 
used  in  estimation  of  that  c.d.f. 

Denote  ordered  sequences  of  event  times  as  follows: 

{t^}  - the  sequence  of  usage  replacement  times, 

{ t j } - the  sequence  of  inspection  times  and  maximum  time  if  any,  and 
{ t^}  - the  sequence  of  ages  of  surviving  engines. 

As  in  part  3,  iik  and  n^  count  usage  removals  and  inspection  removals,  and  j(t^) 
is  the  index  of  the  last  inspection  prior  to  a usage  removal  at  t^.  Conversely, 
i(t!)  is  the  index  of  the  last  usage  removal  prior  to  inspection  at  time  t ^ . 
Similarly,  k(tj|)  and  j(t^)  are  indexes  of  the  last  usage  removal  time  and  inspec- 
tion time  prior  to  the  age  of  a survivor,  t^'.  The  underlying  c.d.f.  of  replace- 
ment times  is  again  assumed  to  have  mass  only  at  observed  removal  times  { t.}  and 
{ t j } , and  these  masses  are  denoted  by  f(").  With  this  notation,  the  likelihood 
function  for  a progressively  censored  sample  from  the  multi-risk  model  (1)  is 
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(usage  removals) 


\+l  j(t.) 

l = c-  n [f (t . ) n (i-p .) ] - 

1=1  1 j=i  3 


k 3-1  i(t  ) "j 

• n [p.  II  (l-p  ) (1-  Z 3 f(t.))]  . (inspection  removals) 

j=l  J i=l  1 i=l  1 


Vi 

•[(1-  z f(t  )) 

1=1 


k 

n 

3=1 


(i-Pj)]" 


k+l 


(max.  time  removals) 


N"  Nk+l-Mk+lr  kK>  3 CO 

k+l  *c+l[ (1-  p f(t  ))  nh  (1_p  )] 

■ 11  k=l  K 3=1  3 

h=l  K 1 31 


(survivors) 


The  log  likelihood  is 


^k+l  j(t  ) 

In  L = In  C + Z [in  f(t.)  + I 1 In  (l-p.)]  + 


i=l 


3=1 


k r i_1  i ( t '. ) 

+ Z n.  [in  p.  + Z ln(l-p. ) + ln(l-  Z 1 f(t.))]  + 

3=1  3 3 i=l  1 1=1 

Vl  k 

+ n [ln(l-  Z f(t.))  + Z In (l-p.)]  + 

. i 1 .1 

1=1  j=l 

N-\+rK 


+ z 

h=l 


k+!  k(t")  j(t") 

[ln(l-  Zh  f(t,  ))  + Zh  In (l-p  . ) ] 


k=l 


j=l 


The  maximum  of  the  log  likelihood  function  will  be  derived  by  calculus  subject 
to  the  usual  constraints  on  c.d.f.s. 

The  first  result  to  be  observed  is  that  the  masses  f(-)  are  all  equal  be- 
tween neighboring  inspection  times  and/or  survivors'  ages  (within  time  intervals 
where  the  numbers  of  survivors  and  inspection  removals  remain  unchanged) . For 
example,  subsequent  to  the  last  inspection  and  last  survivor's  age,  the  derivative 
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is  the  same  for  all  usage  removals  (if  any).  Similarly,  prior  to  any  inspections  and 
survivors,  the  derivative 


' *n  _L  1 
«£(tj  " f(t,> 


k+1  n 

. E _L 

1 = 1 i ( t '. ) 

J 1-  U f(t.) 


N-\+r“k 


k(t") 

kSl^V 


is  the  same  for  all  such  usage  removals.  The  generality  of  the  assertion  is  suggested 
by  the  form  of  the  derivative  subsequent  to  the  last  inspection  but  prior  to  the  last 
survivor’s  age.  (Assume  no  usage  removals  after  the  last  survivor's  age.) 


6 In  L 
6f(t~) 

l 


1 k+1 

f(tr>  ' X 


N-M,  ,-N,  +1 

k+1  k+1 


l-  s f(t.)  l-  r f(t.) 
i=l  1 i=l  1 


The  last  two  terms  combine  to  give  a numerator  which  is  the  number  of  survivors  after 
the  last  usage  removal,  and  the  combination  is  the  same  for  all  such  removals.  In 
order  to  give  the  general  form  of  the  derivatives,  new  notation  must  be  intro- 


duced; k(t!)  is  the  index  of  the  last  usage  removal  prior  to  inspection  at  t ! , h(t^) 

is  the  index  of  the  last  survivor  time  prior  to  the  i*"^1  usage  removal  at  time  t^, 

and  h(t'.)  is  the  index  of  the  last  survivor  time  prior  to  the  j*"*1  inspection  at  time 

t'. . The  general  form  of  the  derivative  is 
1 

N-M,  ,-N,  , 

x r i k+1  lc+1  k+1  , 


6 Jn^L  1 k+1  _nj  k+1  k+1  1 

k=l  k k=l  k 


for  all  usage  removal  i idexes  i between  neighboring  survivors'  ages  and/or  inspection 
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It  is  very  easy  to  show  that  the  estimator  of  the  inspection  time  removal 
probability  has  the  same  form  as  in  section  3.  Set  the  derivative  with  respect  to 


Pj  of  the  log  likelihood  equal  to  zero. 


n.  k+1  n.  N \+l  NP 


= 0 = -J--  S ■= — — - Z 

‘ Pj  Pj  i-j+1  Pj  h=h(t* ) 


The  1-p  terms  factor  out  of  each  sum  leaving 


^ — i - [--  ^ ■]  • [number  of  survivors  after  the  j inspection] 

P i P i 


The  solution  is 


Pj  = n^  / (number  of  survivors  to  the  j inspection) 
as  in  section  3. 

The  relation  (10)  for  the  partial  derivatives  can  be  rewritten  in  the  same 
form  as  equation  (3)  in  the  previous  section,  and  consequently,  it  has  the  same 
form  of  solution.  The  sample  information  must  be  re-expressed  in  order  to  show  the 
relation  between  (10)  and  (3)  and  to  show  the  solution,  the  maximum  likelihood 
estimate  of  F(t) . Combine  the  ordered  sequences  of  inspection  and  survivors' 
ages  into  a single  sequence  { t^},  l = 1 ,2 , . . . . Define  the  new  sequence 
{ r^}  as  the  sequence  of  numbers  of  usage  removal  times  between  successive  inspec- 
tion times  and/or  survivors'  ages. 

It  was  shown  earlier  that  the  estimates  of  f(t..)  are  constant  between 
successive  survivors'  ages  and/or  inspection  times,  and  these  constants  will  be 
denoted  f . The  relation  (10)  can  now  be  rewritten  as 


i 

l 

fi  i-i 


I riA 


and  this  is  the  same  as  equation  (3). 
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In  order  to  express  the  solution,  a generalized  form  of  the  sequence  of 
estimators  { Pj, } must  be  defined  which  includes  the  sequence  of  inspection  time  re- 
moval probabilities.  Let 


number  of  survivors^to  the"~j"fh  inspection  for  l = j inspection 

index,  j=l ,2 , . . . ,k+l 


number  of  survivors  to  the  jth 
survivors'  age 


for  l = j survivors' 
age,  j=l,2,... ,N-Mk+1-K-1 


so  the  maximum  likelihood  estimators  may  be  expressed  as 

f,  = J n [— ],  i-1,2 n-mu+1.  (11) 

Hie  numbers  of  survivors  to  an  inspection  or  to  a usage  removal  can  be  expressed 
as  follows.  The  number  of  survivors  still  in  operation  at  the  time  of  the  in- 
soection  is 


I!  - S n - Em.  - I r , j=l ,2 , . . . ,k+l 
1=1  1 i=l  1 1=1  * 


where  the  index  8 ranges  over  inspection  and  survivor's  ages  up  to  inspection  j or  the  last 
survivor's  age  prior  to  the  jth  inspection  at  time  t j . The  upper  limit  l(t!)  is  the 
index  of  the  jth  inspection  in  the  sequence  of  inspections  and  survivor's  ages.  The 
number  of  survivors  still  in  operation  at  the  time  of  the  jth  usage  removal  may  be 
expressed  as 


j+i  (tj 


N-  Z J(m  « ) - (j-i(t  ))-  s + j-1,2 H. 

i=l  1 1 3 £-1  k 1 


where  i ( t j ) is  the  index  of  the  last  inspection  prior  to  the  j usage  removal. 

Again,  the  index  2.  ranges  over  both  usage  removal  times  and  inspections. 

A function  of  a maximum  likelihood  estimator  is  also  a maximum  likelihood  estimator,  so 
when  r and  p^  are  substituted  into  the  formula  for  the  c.d.f.  F(t)  (1),  the  maximum 
likelihood  estimator  should  be  obtained.  To  define  the  result,  called  the  "product 
limit"  estimator  by  Kaplan  and  Meier  [1],  order  all  usage  removal  times,  inspection 
times,  and  survivors'  ages  into  one  sequence  { t * , r=l,2,...,N.  The  product  limit 
estimator  is 

F(t)  = 1 - n (N-r) / (N-r+1) 
r 

where  r takes  on  values  only  at  usage  removal  times  and  inspection  times  prior  to 
or  at  time  t.  For  example,  if  the  following  realization  had  been  obtained  for  a 
sample 


time 

C1 

J 

c5 

max 

removal 

condition 

Usage 

Survivor 

Usage 

Inspection 

Survivor 

Max. 

Time 

number 

removed 

1 

0 

1 

2 

0 

2 

F(t) 

.125 

.125 

.271 

.562 

.562 

1.0 

The  product  limit  estimator  of  F(t)  in  the  Figure  3 is  obtained. 


Figure  3.  The  Product  Limit  Estimator. 


Figure  4 shows  the  maximum  likelihood  estimator  of  F^(t)  the  usage  removal 
distribution  for  the  same  data.  The  jumps  between  inspection  and  between  sur- 
vivors' ages  will  be  of  constant  height.  The  estimator  of  inspection  removal 
probability  is  (now  carrying  index  l = 2 because  it  is  second  in  the  sequence  of 
survivor  times  and  inspection  times) 

p2=  2/6  = .333. 

Also  needed  are  the  p^  terms  for  the  two  survivor  times  t2  and  t^ 
px  = 1/7  = .143  and  p3  = 1/3  = .333. 

A 

From  (11)  the  jumps  at  usage  removals  are  f = 1/8  = .125  for  the  jump  at  t^ , and 
f2  = (1/8) (1-. 143)”1  = .146  for  the  jump  at  tj. 
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APPENDIX  B 

VARIANCE  REDUCTION  FOR  RENEWAL  AND  REPLACEMENT 
PROCESS  SIMULATION 


SUMMARY 


SUMMARY 


Unless  component  lifetimes  are  exponentially  distributed,  analysis  of 
renewal  processes  requires  transform  methods  or  simulation.  If  components 
are  not  all  identical  or  new,  if  the  process  repeats  for  several  period  of 
time,  or  if  there  are  several  simultaneously  operating  components  supplied 
from  the  same  stock  of  spares,  then  determination  of  replacement  require- 
ments can  only  be  done  by  simulation  or  normal  approximation.  The  use  of 
complementary  antithetic  random  variates  in  the  generation  of  component 
lifetimes  reduces  the  variance  of  the  sample  mean  number  of  replacements 
in  a multiple  component,  multiple  period  replacement  process  simulation 
where  the  components  and  spares  may  be  either  new  or  used. 


MATHEMATICAL  SYMBOLS 


SYMBOLS 


DEFINITION 


PAGE 


t,t.j  time  value  or  age  value 

F(t)  cumulative  distribution  function  of  the 

random  variable  component  lifetime 

R(t)  the  random  variable  residual  lifetime 

given  age  at  installation  is  t time  units 

GI/G/1  the  gueuing  notation  for  a single  server 

queuing  system  with  arbitrary  independent 
arrival  distribution  and  arbitrary  service 
distribution 

N ( t ) the  number  of  replacements  required  to 

obtain  t operating  hours  from  some  identical 
components  operated  sequentially 

l summation 

S(J)  operating  time  required  of  Jth  component 

and  its  replacements 

M the  number  of  simultaneously  operating 

components 

N simulation  run  size 


n 

X,  X. 

Var 

Cov 

U,U. 

x 

e 

In 

N](t) 

{ } 

F'](U) 


set  intersection 

the  random  variable  lifetime  of  a component 

variance 

covariance 

random  numbers  between  0 and  1.0 

lambda,  the  parameter  of  an  exponential 
di stribution 

the  Naperian  base  2.7183 
natural  logarithm 

arbitrary  Lebesgue  inteqrable  functions 

set  containment  symbol 

The  replacement  counting  function  with  a 
stock  of  j spares,  i is  a run  index 

sequences  of  functions  will  be  surrounded 
by  parentheses 

the  inverse  transform  of  F(t) 
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SYMBOL 


DEFINITION 


PAGE 


1 im 

limit 

13 

oo 

infinity 

13 

i 

integral 

13 

E 

mathematical  expectation 

13 

a particular  value  of  a random  number 

13 

1 

used  to  separate  an  event  statement  about 
a random  variable  from  a condition  on  that 
event  in  probabilities  and  expectations 

13 

= 

identity 

14 

Max 

maximum 

14 

Njs(« 

the  replacement  counting  function  with  j 
operating  components  and  initially  s spares, 
k is  a run  index 

16 

Ni(v 

the  counting  function  for  period  j when  t, 
operating  hours  are  required,  i is  a run  J 
index 

18 

M 

mean 

22 

V 

variance 

22 

C 

correlation  coefficient 

22 

VSM 

variance  of  sample  mean 

22 

a 

a parameter  of  the  Wei  bull  distribution 

22 

6 

a parameter  of  the  Weibull  distribution 

22 

H ( J ) 

the  number  of  replacements 

24 

N(I) 

the  cumulative  number  of  replacements  on 
several  runs 

24 

M(I) 

the  cumulative  squared  number  of  replace- 
ments on  several  runs 

24 

ICOV(I) 

the  cumulative  cross  products  of  numbers 
of  replacements  on  several  runs 

24 

CCOEF(I) 

pooled  sample  correlation  coefficient 

25 

VREPL(I) 

pooled  sample  variance  of  the  number  of 
repl acements 

25 

VMEAN(I) 

pooled  sample  variance  of  the  mean  number 
of  replacements 

25 
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VARIANCE  REDUCTION  FOR  RENEWAL  AND 
REPLACEMENT  PROCESS  SIMULATION 


I . INTRODUCTION 

Many  maintenance  systems  call  for  the  replacement  of  operating 
items  at  failure.  It  is  necessary  to  predict  expected  replacement 
requirements  for  reordering  spares  and  for  planning  of  replace- 
ment workload.  It  is  also  useful  for  planning  spares  purchase 
policy  to  know  the  upper  quantiles  of  replacement  requirements  to 
predict  the  probability  of  running  out  of  spares.  The  usual  type 
of  information  available  is  the  cumulative  life  distribution  function 

F(t)  = P[operating  life  <_  t]  t>0  (1  ) 

based  on  accumulated  operating  time  at  failure.  If  a single  operating 
component  is  replaced  at  failure  with  new  identical  components  which  have 
independently  distributed  lifetimes,  then  the  replacement  counting 
variable  N(t),  the  number  of  replacements  required  to  achieve  t operation  time 
units,  is  a "renewal  process,"  and  the  properties  of  Nit)  are,  theoreti- 
cally at  least,  computable,  Cox  [2],  If  a number  of  components  are 
operating  simultaneously,  then  the  number  of  replacements  required 
to  achieve  t operating  time  units  from  each  is  known  as  a superposition 
of  renewal  processes,  and  except  when  the  component  life  distribution 
is  exponential,  only  asymptotic  results  are  known.  In  fact,  the 
only  computationally  convenient  result  for  renewal  processes  and 
their  superposition  is  their  asymptotic  normal  distribution.  These 
circumstances  often  call  for  simulation  of  renewal  processes  for 
prediction  of  replacement  requirements.  Kay  [4]  has  simulated  the 
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upper  quantiles  of  a renewal  process  with  Weibull  life  distribution 
and  compared  the  results  with  the  normal  approximation.  The  normal 
approximation  did  not  appear  to  be  very  accurate  and  tended  to  over- 
estimate actual  replacement  requirements  for  fairly  small  values  of 
operating  hours,  less  than  four  mean  lives.  On  the  other  hand,  the 
simulation  itself  takes  a substantial  amount  of  computer  time  to 
obtain  an  accurate  value  of  the  expected  number  of  replacements, 
especially  for  rare  values  of  the  numbers  of  replacements  such  as 
the  upper  quantiles. 

In  many  maintenance  systems,  the  replacements  are  not  new  but 
used,  or  the  replacements  may  be  a mixture  of  new  and  used  items. 

The  cumulative  distribution  function  of  the  residual  lifetime  R(t) 
given  age  at  installation  t is  the  conditional  distribution 

P[R(t)<_  u]  = [F ( t + u)  -F(t)]  [l-F(t)]-1  u,  t > 0 • 

There  are  no  analytical  results  for  this  more  general  replacement 
process  unless  the  component  life  distribution  is  exponential,  and 
the  memoryless  property  of  the  exponential  distribution  gives  the 
residual  life  the  same  exponential  distribution  as  a new  item, 
regardless  of  age  at  installation.  Simulation  of  this  more  general 
replacement  process  is  necessary  to  obtain  results  about  the 
expected  number  of  replacements  and  its  quantiles. 

The  technique  of  antithetic  variates  was  first  applied  to  simula- 
tion by  Hammersley  and  Handscomb  [5]  for  reducing  the  variance  of 
the  sample  average  output  of  the  simulation  of  the  mean  of  a distri- 
bution. The  technique  has  since  been  applied  to  simulation  of 
61/6/1  service  system  by  Mitchell  [6]  and  to  determinations  of  the 
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critical  path  length  distribution  in  project  networks  with  random  task 
times  by  Burt,  Gaver,  and  Perlas  [1],  but  without  the  order  of  magnitude 
reduction  of  variance  of  the  sample  averaqe  experienced  when  antithetic 
variates  are  applied  to  simulation  of  replacement  processes. 


2 . Simulati on  of  a Replacement  Process 

The  simulation  of  a single  component  replacement  process  follows 
easily  from  the  relation  between  number  of  replacements  W(t)  required 
to  achieve  t operating  hours  and  the  residual  lives  R(t-j)  R^)^,  ... 
for  the  initial  item  and  its  replacements  given  their  ages  at  installa- 
tion t-j , ...  That  relation  is 

P[N(t)=0]  = P[R(t1 ) > t] 

P[N(  t)  = l ] = P[R(t1)<_t  n R(t1)+R(t2)>t] 

etc.  , 

and  in  leneral 

n n+1 

P[N(t)=n]  = P[  I R(t. ) < t n l R(t. ) > t]  (2) 

i=l  1 i=i  1 

To  simulate  the  replacement  process,  generate  the  residual  life  of 
the  first  item,  R(t^).  If  it  exceeds  the  operating  time  t,  count  N(t)=0. 

If  not,  generate  the  residual  life  of  the  second  item  R(t^)  and  count 
N(t)=l  if  R ( t^ )+R ( t^ ) exceeds  t.  if  not,  continue  until  the  sum  of 
residual  lives  exceeds  t.  This  procedure  may  be  repeated  as  often  as 
desired  to  obtain  an  estimate  of  the  mean  number  of  replacements  from  the 
average  of  all  simulations  and  estimates  of  quantitiles  from  sample  auantiles. 

The  simulation  of  a multiple  component  replacement  process  is  simile 
except  that  the  total  operating  time  t is  a vector  of  operating  time  re- 
quirements from  each  component  and  its  replacements.  The  number  of  re- 
placements is  the  sum  of  all  failed  components  and  spares.  This  simulation 
is  shown  in  Figure  1. 


114 


! 


Several  program  details  may  be  supplied  by  the  writer  of  the 
simulation  program.  The  sequence  of  replacements  may  be  specified  or 
a random  replacement  Dolicy  may  select  from  all  available  replacements. 
The  ages  of  replacements  may  be  given  or  generated  from  some  age  distri- 
bution. The  amount  of  operating  time  required  of  each  component  may  be 


specified  as  S(J),  J = 1,  2,  ...,M  for  the  M components.  It  is  no  great 
difficulty  to  include  a resupply  network  that  returns  failed  engines  to 
stock  after  a suitable  repair  delay. 


3 . Application  of  Antithetic  Variates  to  the  Single  Component  Replacement 
Process  Simulation 

The  computer  program  required  to  simulate  the  multiple  component  re- 
placement process  with  high  accuracy  may  strain  the  capacity  of  avail- 
able computers.  One  method  to  reduce  the  run  size  required  to  achieve 
a fixed  level  of  accuracy  is  the  method  of  antithetic  variates,  Hanmersley 
and  Handscomb  [5].  This  method  has  been  used  in  statistics  but  has  re- 
ceived little  use  in  simulation  of  replacement  processes. 


The  measure  of  accuracy  that  may  be  improved  by  the  application  of 
antithetic  variates  is  the  variance  of  a sample  statistic.  For  example, 
the  expected  value  of  a random  variable  may  be  simulated  by  generating 
several  values  of  the  random  variable  X. , i = 1,  2,  ...,  N and  averaging 
them.  The  sample  average  is  an  unbiased,  consistent,  and  otherwise  well 
behaved  estimator  of  the  expected  value.  Its  variance  is 


Var  X 


Var(  l X./N)  = ^ 

i=i  1 r 


I Var  X.  + 2 I I Cov(X  ,X  ) 
‘ i=l  1 i <j  1 J 


(3) 


The  covariance  term  enters  into  the  calculation  if  the  random  variables 

X-j , X2>  ...,  Xf)  are  generated  so  that  they  are  correlated.  The  "trick" 

of  the  antithetic  variates  technique  is  to  generate  them  so  that 

Cov(X.,X.)  <_  0 and  thereby  reduce  the  variance  of  the  sample  mean  below 
N 1 3 2 

l Var  X./N  without  biassing  the  sample  mean. 
i=l  1 

Two  methods  exist  for  creating  negative  correlation  among  random 
numbers  used  to  generate  random  variables.  They  are  to  use  the  complements. 
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(1 -random  number),  and  to  generate  half  the  required  set  of  random  num- 
bers and  use  them  in  reverse  order  for  the  second  half. 

Both  methods  were  tested,  and  the  former  seems  more  promising 
for  simulation  of  renewal  and  replacement  processes.  The  computer  pro- 
gram is  shown  in  Figure  2.  One  reason  for  this  preference  is  that  many 
random  numbers  must  be  generated  to  give  one  value  of  the  number  of 
removals  according  to  the  relations  (2).  Reversing  the  sequence  of 
random  numbers  used  to  generate  the  R(t^)  in  the  next  run  may  give  inde- 
pendent values  of  N(t)  if  the  number  of  replacements  is  small  relative 
to  the  number  of  random  numbers  generated.  For  example,  suppose  a test 
program  generates  50  random  numbers  for  each  run,  U(l),  U(2),  ...,  U(50), 
which  are  used  to  generate  50  residual  lives.  If  t is  small,  then  the 
number  of  replacements  may  also  be  small,  say  5.  On  the  second  run, 
the  sequence  of  random  numbers  will  be  used  in  reverse  order,  U(50), 
11(49),  ...,  U(l),  and  N(t)  will  probably  involve  only  the  first  few 
values  of  the  reversed  sequence,  independent  of  U(l),  U(2),  ...,  U(5). 

The  variance  reduction  property  will  be  proved  for  the  single  com- 
ponent replacement  process  and  for  several  methods  of  aenerating  the 
residual  life  times  of  the  components,  the  inverse  transformation  method 
and  a method  of  acceptance  for  generating  conditional  lifetimes.  The 
variance  reduction  property  for  mul ti -component  replacement  processes 
will  be  proved  in  the  next  section. 

The  inverse  transformation  method  is  based  on  the  fact  that  for  any 
strictly  increasing  cumulative  distribution  function  F(x),  the  random 
variable  F(X)  has  the  uniform  distribution  on  the  interval  from  zero  to 
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one.  The  inverse  transform  of  a uniform  random  number  U is  denoted 
F ^ ( U)  and  is  a random  variable  with  cumulative  distribution  function 
F(x). 

This  procedure  for  generating  random  variables  can  be  carried 
out  on  a graph  of  F(x)  by  reading  across  from  a random  number  on  the 
vertical  axis  to  the  curve  F(x)  and  down  to  the  horizontal  axis  for 

a value  of  the  random  variable.  Figure  3.  This  method  may  be  adapted 
for  generating  random  variables  from  a discrete  distribution. 


FIGURE  3 ; INVERSE  TRANSFORMATION  METHOD 
OF  GENERATING  RANDOM  VARIABLES 

The  method  of  inverse  transformation  is  applicable  only  for 
families  of  cumulative  distributions  which  can  be  solved  exDlicitly  for 
the  inverse  F~^(U)  (or  for  distributions  entered  in  tabular  form  and 
random  variables  generated  as  in  Figure  3).  For  example,  the  ex- 
ponential cumulative  distribution  function  F(x)=l-exp(-Ax)  may  be 
inverted  to  generate  an  exponential  random  variable  X as  a function  of 
a uniformly  distributed  random  number  U,  X=-ln(l-U)/x.  This  method  may 
also  be  used  to  generate  Weibull,  Cauchy,  logistic,  Rayleigh,  and  of 
course,  uniform  random  variables. 
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The  variance  reduction  property  to  be  proved  is  : 

Theorem  1 : The  use  of  complementary  antithetic  random  variates  reduces 

the  variance  of  the  sample  mean  number  of  replacements  in  simulation 
of  a single  component  replacement  process  when  component  lifetimes  are 
generated  by  the  inverse  transformation  method. 

Proof:  According  to  (3),  the  variance  can  be  reduced  by  generating 

correlated  replacement  counting  functions  so  that  for  every  pair  of 
runs,  the  covariances  of  the  numbers  of  replacements  CovlN1 (t) ,NJ ( t) ) 
are  non-positive  thereby  reducing  the  variance  of  the  sample  mean  from 
N runs  to  less  than  or  equal  to 

? N 

(1/N  ) l Var  N 1 ( t ) 
i=l 

obtained  by  generating  the  replacements  in  each  run  from  independent 
random  numbers.  It  will  later  be  shown  that  for  at  least  one  pair  of 
runs,  the  covariance  is  strictly  negative. 

The  proof  of  non-positive  covariance  is  to  apply  the  following 
result  from  Mitchell  [6]  to  replacement  processes  with  a limited  number 
of  spares:  "Let  f and  g be  real  valued  Lebesgue  integrable  functions  of 

k 

XeR  . Suppose  for  each  argument  X.,  l^j_fk,  either  f is  non-increasing  and 

3 

g is  non-decreasing  or  f is  non-decreasing  and  g is  non-increasing  in  Xj. 
Let  U= (U, ,LL, . . . ,1^)  be  uniform  on  [0,11.  Then  Cov(g(U) ,f (U) )<0. " 
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Define  the  sequences  { N 1 ( t ) } and  {N?( t ) } for  j = 1,  2, 

J J 


follows : 


N-J  ( t)  = 0 

, ( 0 if  F_1  (U, ) > t 

No(t)  = ' 

* ( 1 if  F 1 (Un ) i t 


0 if  F-^ ( U-j ) > t 

i]3(t)  = 1 if  F-1  (U-| ) 1 t n F-1  (U-j ) + F_1(U2)  1 t 
2 if  F- 1 ( U-j ) + F-1(U2)  i t 


and  N^(t)  = 0 


N|(t)  = 


0  if  F-  ( 1 — U] ) > t 


1  if  F"1  (1-^)  l t 

0 if  F-1  (1 -U-, ) i t 

1 if  F-1  (1  -U-j ) £tn  F-1  (1-U] ) + F_1(l-U2)  > t 

2 if  F_  1 ( 1 - U-j ) + F_1(l-U2)  < t 


The  Nj(t)  are  counting  functions  which  tell  the  number  of  replacements 

required  given  a supply  of  j-1  replacement  items.  The  Nj(t)  sequence 

is  generated  by  inverse  transformation  F ^(11)  from  ranoom  numbers 

U-, , U , U.  , and  the  N.(t)  sequence  is  generated  from  the  comple- 
^ J ^ J 

ments  1-U-|,  1-U2>  ...»  1-lL 

First  it  will  be  shown  that 


Cov(N](t),N2(t))  < 0 
J J 


■ 


for  all  j and  then 

1 im  Cov(n!(  t)  ,N2(  t) ) = Cov(N^ ( t) ,N2(t) ) ^0 
j-*~  J J 


follows  directly.  To  apply  Mitchell's  result,  it  must  be  shown  that 

1 2 
N j ( t ) is  nondecreasing  in  each  argument  , U^,  ^ and  N ^ ( t ) 

is  nonincreasing  in  the  same  arguments.  This  occurs  because  a large 

random  number  generates  a large  random  1 i f e t i me  F ^(U^)  resulting 

in  a small  N.(t)  and  conversely  for  N.(t)  because  large  U,  results  in 
3 J ■ 

2 

small  1 -Ui and  small  F (1-U^)  and  a larger  number  of  replacements  Nj(t). 
The  fact  that  all  the  covariances  are  bounded  allows  the  result 


lim  Cov(n1 ( t) ,N2(t) ) = Cov(lim  N^(t) ,lim  N2(t)) 
j-*»  J J j^co  J j-K»  J 

= CovCN1 (t),N2(t)) 

the  covariance  of  replacement  processes  with  unlimited  supplies  of 
replacements.  Since  every  term  in  the  sequence  of  covariances  is  non- 
positive, so  is  the  limit 

Cov(N1(t),N2(t))  < 0 . 


In  order  to  show  there  will  be  some  improvement  when  complementary 
antithetic  random  variates  are  used,  we  express  the  covariance  condi- 
tioned on  the  value  of  the  first  random  number  U-j  which  has  the  uniform 
density  on  the  unit  interval. 

Cov(N1(t),N2(t))=  /'  E [ (N1 (t)  - EN(t))(N2(t)  - EN(t))|  U,  = u]du 

1 2 

where  EN(t)  = EN  { t)  = EN  (t).  Assume  there  exists  a random  number  u 
close  enough  to  1.0  so  that  F_1(u)  is  larger  than  t.  Then  N^t)  0 
when  U-j  > u and  the  covariance  integral  is 


« 


f 
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« 

1 • 

Cov(N1  ( t)  ,N2(  t)  )=  / -EN(t)E[N2(t)-EN(t)|U,«  u]du+ 
u 1 

+ J~  E[N]  ( t)-EN(t) ) (N^(t)-EN(t) ) [ U,=u]du. 

0 1 

The  first  term  will  be  strictly  negative  because  U^u  means  1-U-jis 
small  causing  E[N2(t)|U1=u  ] to  exceed  EN (t) . The  second  term  will 
remain  non-positive  because  N (t)  is  larger  than  usual  and  N (t) 
smaller  under  the  condition  U-j  = u <u.  Thus  some  reduction  in  the 
variance  of  the  sample  mean  can  be  expected  when  complementary  anti- 
thetic variates  are  used  on  alternate  simulation  runs. 

This  first  theorem  is  contingent  on  generating  the  component  lives 
by  inverse  transformation.  When  comoonents  have  already  accumulated 
operating  time  prior  to  installation,  their  residual  operating  times 
must  be  generated  from  the  conditional  life  distribution 
P[R(t)<_U|Life>t.]  • 

The  inverse  transformation  may  again  be  used  and  the  conditional  inverse 
must  be  recomputed  for  each  different  age  ti . The  variance  reduction 
property  will  still  hold  when  antithetic  random  variates  are  generated 
by  inverse  transformation  from  the  conditional  life  distribution  above. 

The  same  method  of  proof  yields: 

Corollary  1 : The  use  of  complementary  antithetic  random  variates  reduces 

the  variance  of  the  sample  mean  number  of  replacements  in  simulation  of  a 
replacement  process  when  the  component  residual  lifetimes  are  generated  by 
the  inverse  transformation  method  from  their  conditional  life  distributions 
given  their  ages  at  installation. 


f 
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4 . Multiple  Component  Replacement  Processes 


Often  the  replacement  process  will  include  more  than  one  operating 
item  and  a stock  of  replacements  for  any  item  as  it  fails.  If  all  re- 
placements are  new,  then  the  number  of  replacements  is  obtained  from 
superDOsi tion  of  renewal  processes  for  each  ooerating  item  and  some  in- 
formation about  the  distribution  of  the  number  replacements  is  known, 

Cox,  chapter  6 [2],  If  the  replacement  items  and  perhaos  the  operating 
items  are  not  new,  simulation  is  the  only  way  to  estimate  the  distribu- 
tion of  the  number  of  replacements  and  other  properties.  The  simulation 
program  is  shown  in  Figure  1.  It  can  be  adapted  to  simulate  any  policy 
that  selects  the  sequence  of  installation  for  spares  including  a 
random  replacement  policy. 

Complementary  antithetic  variates  can  again  be  used  to  reduce  the 
variance  of  the  sample  mean  number  of  replacements  for  the  mul ti -component 
replacement  process.  The  application  is  to  simulate  the  replacement 
process  to  obtain  the  replacement  requirements  to  achieve  the  desired 
operating  time  t for  all  items  and  replacements.  Save  the  random  numbers 
used  in  generating  item  life  times  and  use  their  complements  in  the  next 
simulation  run  to  generate  component  lives  in  the  same  order  as  they  are 
used  in  the  first  run.  The  variance  reduction  property  will  be  proved 
under  the  assumption  that  the  item  lives  are  generated  by  the  inverse 
transformation  method; 

• Theorem  2:  The  use  of  complementary  antithetic  random  variates  reduces 

the  variance  of  the  sample  mean  number  of  replacements  in  simulation  of  a 
multiple  component  replacement  process  when  the  item  lives  are  generated 
by  the  inverse  transformation  method. 

I 

* 

S 
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Proof:  The  method  is  as  in  Theorem  1 to  limit  the  number  of  replacements 

available  and  use  Mitchell's  Theorem  to  prove  the  non-positivity  of  co- 
variance  between  two  antithetic  counting  functions.  Then  the  limit  on 
replacements  can  be  removed. 

Define  the  item  lifetimes  for  a replacement  process  that  starts  with 
j simultaneously  ooeratinq  items  and  k spares  as  , and  the  antithetic 

sequence  of  lives  as  X2i  , i = 1,  2,  . ..,  j+s.  The  counting  functions 
obtained  from  the  two  seauences  will  be  obtained  as 


0 if  X^.  > t i=l , 2,  . . . , j 

II  if  only  one  X^.  £_  t,  i =1 , 2,  . . . , j 

and  X. survives  past  time  t 
kl+1 

j+s if  all  s replacements  are  needed, 
k=  1,2.  Since  the  inverse  transformation  method  is  used  to  generate 
all  random  lifetimes,  the  counting  functions  can  be  expressed  as  functions 
of  j+s  random  numbers, 

n]sU)  • f(u, , u2,  ....  uJ+s) 

and 

Njs(t)  = gO-U-j , 1-U?,  ...,  l-uj  + s). 

Mitchell's  Theorem  can  again  be  used  to  show  that 
Cov(Njs  (t) , N jS  (t) ) <_  0 

and  consequently  in  the  limit  as  s the  number  of  spares  increases.  The 
proof  of  Theorem  1 can  be  adaDted  to  show  these  will  be  strictly  negative. 


As  in  the  previous  oroofs,  the  effect  of  a decrease  in  will  be 


decreases  to  less  than  t,  requiring  a replacement.  The  opposite  effect 

occurs  on  Njs(t)  because  is  generated  from  1-1^  and  is  larger  if 

1 2 

U-|  decreases.  Then  Njs(t)  and  li.s(t)  move  in  opposite  directions  in 
, response  to  changes  in  or  for  that  matter  any  IK  i = l , 2,  ...,  j+s 

although  the  effect  may  not  be  as  large  for  large  i because  the  cor- 
responding spare  item  may  never  be  needed. 

The  remainder  of  the  proof  is  to  let  the  number  of  spares  increase 
to  <», 

tim  Cov  (n!  (t),  (t) ) = Cov  ( N 1 ( t ) , N? ( t ) ) <_  0. 

J J J 

for  the  counting  functions  of  a j item  replacement  process  where  all  j 
items  or  their  replacements  must  operate  for  t time  units  and  the  number 
of  spares  is  unlimited. 

The  complimentary  antithetic  variates  improvement  in  variance  of  the 
sample  mean  number  of  replacements  can  be  extended  to  the  case  where  the 
required  operating  time  for  each  item  is  different,  t is  a vector  of  j 
operating  times  required  from  j items  or  their  replacements.  The  proof 
includes  the  situation  where  not  all  items  and  spares  are  new,as  lonn  as 
the  inverse  transformation  is  used  to  generate  residual  lives. 
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5.  Multiply  Period  Real acenient  Processes 

Often  the  replacement  process  will  proceed  for  many  periods  and 
estimates  of  the  replacement  requirements  in  each  period  are  needed. 

Because  items  which  have  been  operated  in  one  period  are  survived  have 
accumulated  operating  time  when  put  into  service  in  the  next  period,  the 
situation  in  every  period  but  the  first  is  as  in  corollary  1 of  theorem 
1,  the  items  to  be  used  may  already  have  accumulated  a random  age  depend- 
ing on  their  initial  age  if  any  and  their  accumulated  operating  time  in 
earlier  periods.  Corollary  1 of  theorem  1 may  be  aoplied  to  prove  negative 
covariance  is  obtained  if  antithetic  random  variates  are  used  to  generate 
residual  random  lives  for  all  components  given  their  age  at  the  beginning 
of  each  period  where  they  may  be  used.  This  will  require  regeneration  of 
the  residual  life  of  any  component  used  in  more  than  one  period  from  its 
conditional  life  distribution  given  its  age  at  installation  in  a subsequent 
period. 

For  example,  let  N^t^)  and  N1^),  i = 1 , 2 be  the  numbers  of  replace- 
ments in  periods  one  and  two  requiring  operating  times  t-j  and  t^  when 
generated  from  independent  random  numbers  and  their  complements.  Random 


numbers  U-|  > , . . . , U, 


Hl(t  are  used  to  generate  component  lives 


X11  ’ •••’  Xl,N1(t1)+l  Where 


,1 


(tJ+1 

Component  number  N^'t^+l  survived  to  be  used  in  the  second  period  with 
accumulated  age 
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N ( tj ) 

xn  a 

i=i 

Random  number  U^l ^ ^ is  used  to  qenerate  the  residual  operating  life  in 

period  2 from  an  inverse  transformation  on  the  conditional  life  distribution 

F(a+t)-F(a) 

l^FTal 

and  all  subsequent  replacements  will  be  generated  as  usual.  This  means  that 
the  replacement  processes  in  several  periods  given  the  ages  of  components 
that  may  have  survived  Dartial  use  in  earlier  Deriods  are  conditionally  in- 
dependent because  the  random  numbers  used  to  generate  each  period's  replace- 
ment process  are  different  and  independent. 

The  only  remaining  problem  is  to  insure  that  the  antithetic  replacement 

1 2 

processes  are  also  conditionally  independent.  Let  N=max  (N  (t^),N  (t^))+l 
and  use  random  numbers  U^,  ...  to  qenerate  the  reolacement  process 

and  its  antithetic  process  for  period  2 independently. 

The  results  of  this  extension  can  be  adapted  for  reduction  of  the 
variance  of  the  sample  mean  for  a multiple  component,  multiple  period 
replacement  process. 


6.  EMPIRICAL  TESTS  OF  ANTITHETIC  VARIATES 

In  this  section  results  on  several  types  of  replacement  process  simulations 
with  complementary  antithetic  variates  are  reported.  Although  run  sizes  are  not 
large  and  results  are  somewhat  dependent  on  the  seed  chosen,  the  reduction  in  the 
variance  of  the  sample  mean  number  of  replacements  is  unexpected.  The  proofs  in 
previous  sections  only  show  there  will  be  a reduction  but  do  not  suggest  the 
magnitude  of  the  reduction  which  in  several  cases  was  one  tenth.1  (Mitchell  {6} 
obtained  about  a 20%  reduction.) 

Runs  were  made  for  several  different  combinations  of  new  and  used  components, 
single  and  double  operating  components,  and  several  distributions  of  ages  and 
lives.  In  every  run,  the  unconditional  expected  life  was  0.5,  and  each  component 
and  its  replacements  were  to  operate  for  either  five  or  ten  hours.  Ages  for  used 
components  were  generated  from  the  same  distribution  as  the  component  lives,  but  the 
choice  of  age  distribution  is  arbitrary.  The  total  number  of  components  used 
was  recorded  for  each  run.  The  simulation  program  would  first  make  a number  of 
runs  (20  or  100)  with  independently  generated  random  numbers.  The  results  from 
these  runs  are  listed  in  column  "Base".  Then  an  equal  number  of  runs  were  made 
using  the  complementary  antithetic  random  numbers.  A third  set  of  runs  was  generated 
from  another  set  of  independent  random  numbers  and  the  result  of  these  runs  are 
tabulated  in  the  column  "Independent".  The  means  and  variances  of  the  number  of 
replacements  in  the  Base,  Complementary  and  Independent  runs  are  tabulated.  Then 
the  correlation  coefficient  between  Base  and  Complementary  and  between  Base  and 
Independent  are  listed  in  columns  Complementary  and  Independent.  In  nearly  every 
simulation,  the  correlation  between  Base  and  Complementary  numbers  of  replacements 
is  less  than  -.5  and  the  correlation  between  Base  and  Independent  is  moderate.  The 
last  entry  is  the  variance  of  the  sample  means  computed  from  the  formula  (3)  after 


combining  variances  and  means  in  pairs.  Base  with  Complement  and  Base  with  In- 
dependent. The  simulation  results  are  in  Table  1. 

The  statistical  subroutine  which  was  used  may  be  adapted  to  any  test  of 
effectiveness  of  any  variance  reduction  technique  applied  to  simulation  of  any 
output  statistic.  It  is  flowcharted  in  Figure  4 and  listed  in  FORTRAN  G in 
Figure  5. 


TABLE  1 SIMULATIONS  OF  REPLACEMENT  PROCESSES 
(continued  from  previous  page) 


Input  number  of  runs  K,  operat- 
ing  components,  operating  time, 
and  age  and  life  c.d.f.s. 


Generate  ages  for  all  components 
and  replacement  items  for  use  in 
every  run. 


Compute  cumulative  sums  and  sums 
of  squares  ~ 
N(I)=N(I)+H(J),  M(I)+M(I)+H(Jr 


Compute  cumulative  cross  products 
for  0=1  and  2 and  J=2  and  3 pairs. 
I CO V ( I ) = I CO V ( I )+H( J-l  )*H(  J) 


Compute  means,  variances,  covar- 
iances, and  pooled  variances  of 
sample  means  from  M,  N,  and  ICOV. 

, I . 

/ output  J 
( Stop  ^ 


figure  4,  SIMULATION  TEST  PROGRAM 
134 
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H( 1 ) , H(2),  AND  H( 3)  ARE  THE  NUMBERS  OF  REPLACEMENTS 
REQUIRED  ON  A RUN  USING  BASE,  INDEPENDENT,  AND  COMPLE- 
MENTARY RANDOM  NUMBERS.  N(I)  AND  M(I),  1 = 1 ,2,3,  KEEP 
RUNNING  SUMS  AND  SUMS  OF  SQUARES,  AND  ICOV(I),  1=1,2, 
COLLECTS  CROSS  PRODUCT  SUMS. 

DO  100  1=1,3 
N(  I )=N( I )+H( I ) 

M( I )=M( I )+H( I ) *H( I ) 

100  CONTINUE 

DO  101  1=2,3 

I CO V ( I ) = I COV ( I )+H( 1-1 )*H( I ) 

101  CONTINUE 

VAR(I),  1=1, 2, 3,  IS  THE  SAMPLE  VARIANCE  OF  THE  NUMBER 
OF  REPLACEMENTS.  COV( I ) , CCOE F( I ) , VREPL(I),  AND  VMEAN(I), 
1=2,3,  ARE  POOLED  COVARIANCES,  CORRELATION  COEFFICIENTS, 
SAMPLE  VARIANCES,  AND  VARIANCES  OF  SAMPLE  MEANS  FOR 
BASE/INDEPENDENT  AND  INDEPENDENT/COMPLEMENTARY  PAIRS. 
VARIABLE  "K"  IS  THE  NUMBER  OF  RUNS  IN  EACH  OF  BASE, 
COMPLEMENT,  AND  INDEPENDENT  SIMULATIONS. 

DO  102  1 = 1 ,3 

102  VAR( I )= (M( I )-(N( I )*N( I ) )/K)/ (K-l ) 

DO  103  1=2,3 

COV( I )=( I COV( I )-K*MEAN( 1-1 )*MEAN( I ) )/ ( K-l ) 

CCOEF( I )=COV( I )/SQRT( VAR( 1-1 ) *VAR( I ) ) 

CME AN ( I ) = ( ME AN  f I - 1 ) +MEAN ( I ) ) /2 . 0 

VREPL(I )= (M( I - 1 )+M( I)-( (2*K)*CMEAN(I)**2) )/(2*K-l ) 

103  VMEAN ( I ) = ( VAR ( I - 1 ) +VAR ( I ) + ( 2 *COV ( I ) ) ) / ( 4*K) 

STOP 


Figure  5 STATISTICAL  SUBROUTINE  LISTING 
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